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I. INTRODUCTION 


A composite, by definition, is any two or more materials combined on a 
macroscopic scale to form a useful material [Ref 1]. Today, in one form or 
another, composites are being used on every level of our society. This thesis 
will focus on a class of composites known as "fibrous composites,’ today these 
are typically used in high tech applications, often high performance military 
aircraft. Fiber composites consist of very thin fibers with high aspect ratios and 
high strengths imbedded in a matrix, commonly a polymeric matrix such as 
epoxy. Glass/epoxy is an example of a frequently utilized composite mixture 
that has been in use for many years. 

One of the many advantages this class of material provides is high 
strength-to-weight, though due to a lack of thorough understanding in failure 
mechanics, composites have been limited to secondary load carrying structural 
designs such as skins and control surfaces. In recent years (last 20), composites 


are expanding into a much wider role. 


A. REQUIREMENTS IN ADVANCED STRUCTURES 

Since the introduction of glass/epoxy, advances in production techniques and 
developments in the textile industry have resulted in fiber materials with much 
improved strengths and stiffness. Two of the early break-throughs were 
graphite and boron fibers which have typical strengths and _ stiffness much 


higher than convention ductile isotropic metals. These fibers allow for new 


design limits but also give mse to the requirement of new design technique 
utilizing a maternal of an unconventional material redundancy with anisotropy. 

Employing the fiber/matrix as unidirectional lamina to form multidirectional 
laminate gave the designer the tools to make structural designs subjected to 
multiple loading conditions a manageable task. This idea shifted the attention 
of the designer from the microscopic level to the macroscopic level. 

Utilizing the high strength and stiffness, fibers application for structural 
enhancement were incorporated in many of the military designs in the early 
70’s such as the F-111 aft fuselage structure. Testing of these designs were 
accomplished by component. This method is time consuming and expensive but 
must be done. To gain any type of reliability assurance many tests are 
required. It soon becomes apparent that testing costs could out pace production 


costs. 


B. FUTURE DESIGN REQUIREMENTS 

As weaponry becomes more and more sophisticated the demands on 
composite materials for a wide variety of uses are ever increasing. New designs 
often are only limited by the structural limitations of the materials in use. 
Such restrictions might include the thickness of a wing or stabilizer, the sweep 
of the wing, or the "g" loading. Lately new design specifications have included 
the use of composites for the radiation energy absorbing properties. Though the 
Navy’s interest in composites spans a multitude of uses, one special interest is 
in the area of structural enhancement, such as rocket motor cases, pressure 


vessels for submarine flasks and jet aircraft pilot ejection seat. 


In the last few decades composites have advanced dramatically, and are 
beginning to provide improved alternatives to more conventional methods of 
design normally utilizing isotropic alloys. Fiber reinforced plastics have 
attracted a large amount of attention. Improved techniques in the textile 
industry have lead to more uniform fibers from bases of precursors for high 
performance graphite’s with strengths and stiffness exceeding high yield steels. 

As utilization goes up the increased use of new fibers in composites and the 


cost of system failure testing demand better methods for reliability evaluation. 


II. BACKGROUND 


A. COMPOSITE LOAD SHARING 

A brittle material generally has a higher strength to density ratio but a 
lack of ductility limits the structural uses. Ductility resists crack propagation 
from a region of damage. Metals have a high ductility and do provide practical 
strength and stiffness properties for many structural applications. Fibers do not 
use ductility to resist failure [Ref. 2:pp. 1-19], because ductility is not required 
the fiber can be synthesized from brittle materials. 

The maximum theoretical elastic strain for carbon hydrogen bonds is 
approximately ¢,,, = 0.1, which is one order of magnitude greater than 
present day strong elastic reinforcing fibers. 

Presently fibers can not obtain strengths any where near the theoretical 
ideal strength of the material but are an order of magnitude greater than 
conventional bulk materials. Fiber strength can also be attributed to the small 
fiber diameter, which limits the flaw sizes. 

Composites are geometrically advantageous to homogeneous material such as 
an all metal or all ceramic material, for its ability to utilize the high strengths 
and stiffness of brittle materials while maintaining strength redundancy through 
the matrix. 

The strength interaction of the fiber and matrix is the key to reliability 
effectiveness of a composite. Consider several fibers in parallel without matrix. 
The load is equally shared by each fiber. If one fiber should fail the failed fiber 


is completely ineffective and no longer supports any load. The load of the 


remaining fibers increases equally. With a matrix present (Figure 1) a crucial 
phenomenon occurs. At the point of fracture the more ductile matrix absorbs 
the load through shear stresses along adjacent fibers and back to the broken 
Gber itself. Near the fracture site high shear stress develope in the matnx 
along the broken fiber but quickly 
dissipate as the normal stress 
increases. On the two surrounding 
fibers a noticeable stress increase is 
observed adjacent to the fracture. 
The short distance from the broken 
end of the fiber to the point that 
fiber is carrying a full stress is called 
the ineffective length [Ref. 3] 

The ability of the composite to 


adjust to multiple fiber failures is it's 





strength redundancy, and is Figure 1. Local load sharing in a 
Fiber/Matrix composite. 

dependent on the characteristics of 

both the matrix and fiber. There are several fracture mechanisms all of which 

are a function of the bond strength (shear strength of interface), the matrix 

shear strength and the distribution of the fiber strength. The most critical is 

the fiber strength distribution. The distribution of the fiber determines the 


composite reliability. 


B. STATISTICS AND RELIABILITY 

The complexities of composites make the task of analysis virtually 
impossible. New fibers are being developed every year and attempting to 
generate thorough material testing studies as has been done with conventional 
matenials is not realistic. But composite reliability can be estimated through 
probability studies. Probability of the composite must be based or inferred from 
accurate Statistics of the fiber distribution. 

To determine the reliability of the composite, the strength distribution of 
the fiber must be known. The distmbution can be quantified by strength 
statistics (a Empirical Probability Distnbution Function (EPDF), or Empirical 
Cumulative Distribution Function (ECDF)). For computational efficiency the 
empirical data can then be modeled with an analytical function and thereby 
characterize the distribution by parameters of the function. To define the 
parameters the distribution function must be identified by model type (..e., 
Normal, Weibull, etc.). Using the model Probability plot; the closeness of fit of 
the data can be judge and therefore the appropnateness of the model. Selecting 
the model allows for the calculation of the model parameters. 

The relationship between the fiber statistics and the composite probability is 
represented in a linearized Probability plot (Figure 2). 

Reliability 1s mostly concerned with the lower tail of the distribution, this 
also being the most difficult to define. Because of the redundancy provided by 
the matnx induced load sharing the composite is considerably more reliable 
than it’s constituent fibers in the lower tail. As the fiber statistics fluctuate in 


the extreme lower tail the composite reliability fluctuates with it. This can be 


Probability Plot 
( representation } 


Composite Ul 





Figure 2. Comparison of Fiber and Composite in a 

Probability Plot.. 
shown in a representation of a probability plot in Figure 3. There is a critical 
point where the fiber plot stabilizes, the linear region in this plot determines 
the accurate statistics of the lower bound reliability (worse case) of the 
composite. 

Fibers have a characteristic which helps to identify the analytical 
distribution model that will closely define the empirical statistic. Traditionally, 
the weakest link model [Ref. 4] depicted the fiber as being made of many small 
segments linked together much like a chain. The segments each have an 
intrinsic strength, none being exactly equal. When the fiber is stressed the 


weakest will fail and the total fiber fails. The Weibull distmbution mode! 


defines this weakest link behavior and is commonly used in modeling fiber 


strength. 


1. Weibull Density Distribution Function 


The Two parameter Weibull cumulative distribution function (CDF) 1s 


defined by the following equation. 


F(x) = 1 - exp {-(x/B)*}, x, 8, a >0 (2.1) 
where: x = random variable (stress, load) 
8B = scale 
a = shape 


The parameter alpha defines the shape (spread and skew) of the 


Probability Plot 
( representation } 





Figure 3. Extreme lower tail. 


distribution. Alpha less than 3.5 the 
plot is skewed positively, greater than 
3.5 it is skewed negatively and 
neutral at 3.5. The alpha parameter 
determines the strength scatter of the 
fiber, an increasing alpha defines less 
scatter and higher strength 
uniformity. The beta parameter 
defines the central tendency of the 
distribution and is relatable to the 
mean. Alpha for a fiber type is 


constant for the weak link physics, 


whereas beta changes with the fiber length. An increase in length causes a 


decrease in beta (a weaker fiber). 


The composite reliability can be inferred from the distribution of a fiber 
strength. To gain high reliability with high confidence levels a proportionate 
amount of statistics must be obtained (i.e. for a reliability level of 10° requires 
10° statistics). 

Presently fiber statistics are gathered by single filament testing (Figure 
4). First a single fiber specimen is prepared, since a fiber has such a small 
diameter, (AS-4 diameter ~ 10 micron) it is very fragile. The fiber is placed in 
a cardboard frame. Once the specimen is loaded in the test apparatus the 
cardboard frame is severed by a 


heated wire followed by _ load 


Single fiber testing 


application and interpretation of the 


load deformation data. This is a PEN afi wievsvnll 
Case 


time consuming effort which aan 


produces only one failure’ strength 





statistic. To gain the number of & 
Figure 4. Single filament test 
statistics necessary would take specimen. 
thousands of tests, as an example man-safe reliability is 10° therefore to 
achieve the necessary statistics would require 10’ tests. 
A new approach being explored is to test many fibers in parallel (a 
bundle) and recover the fiber characteristics from the failure test. This would 
provide a large statistical base required and could be accumulated with a much 


higher efficiency. With a validated method the analysis of fiber reliability could 


keep pace with the development of new fibers with reduced effort. 


2. Size Effect. 
The extreme lower tail is difficult to determine accurately for the fiber. 
For a Weibull failure physics and the resulting strength distribution the lower 
tail can be estimated by assuming weakest link [Ref. 4]. The parameter beta 


is a function of the fiber length and is defined by equation 2.2 
B, = 8,0,A,)°° (Zs 


where beta for 1, 
beta for 1, 
Q, 


B, 
B, 
O =e 


A detailed development can be seen in Appendix A. Using this 
relationship parameters for beta can be retrieved analytically where practical 
testing is impossible and the extreme lower tail is estimated. 

This thesis will attempt to measure the lower tail via bundle testing, as 
opposed to the aforementioned estimation, (fiber statistics via single filament 
testing). This investigation explores the methodology needed in testing a bundle 
of fibers to failure. By recording the load and displacement it is expected that 
the underlying distribution of the fibers in that bundle may be extracted. If, 
for example, a bundle has 3k of fibers, the distribution is based on 3k of 
statistics. As the statistical base approaches the complement of the reliability, 


the accuracy of the estimation increases. 
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I. TECHNICAL ISSUES AND RESOLUTIONS 


A bundle of fibers tested to failure presents several experimental obstacles 
which effects the analysis. There are three experimental consequences to 
testing a bundle in parallel that are not experienced in testing a single fiber in 
series. These phenomenon which may disturb the failure distribution, are 1) a 
protective/adhesive enhancement coating called sizing which is placed on the 
filament spools during production, 2) friction which occurs during testing after a 
portion of the bundle has failed and 3) slack which results from specimen 


fabrication. Each of these is discussed in the following paragraphs. 


A. SIZING 

Fibers are packs on spools of bundles, in order to prevent the fibers 
from entangling a small amount of a chemical liquid called sizing is coated on 
the bundle. The effect of this sizing on bundle testing is unknown and must be 
address at sometime. The testing in this experiment was done leaving sizing on 
the bundle. A method to remove the sizing chemical without damage to the 


fibers needs to be studied but is not addressed here. 


B. FRICTION 

During the bundle test prior to fiber failure there is negligible friction, this 
is due to the equal strain between fibers. As the fibers fail the strain of the 
failed filaments return to zero leaving the loose ends entangled in the 
remaining fibers, this places added strain on those fibers. The resulting effect 


shows up as reduced failure strengths and clustering of the these failures. In 
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the relatively short gauge lengths friction is expected to be small and very 
pronounced in the long gauge lengths. The effects of lubricating the bundle 


with oil during testing will be observed to see if the effect is reduced. 


C. SLACK 

During specimen preparation and placement in the test apparatus the 
alignment of the fibers becomes misaligned with unequal length. The effective 
result is an initial nonlinear load/displacement which also effect the failure 
distribution. Another concern is the deformation in the actual load train of the 
test apparatus. This is referred to as system compliance. Compliance was 


found to be a function of the load and can be removed from the resulting data. 


MW 


IV. EXPERIMENTATION 


The experimental setup consisted of a Material Tester (INSTRON 4206) 
and an integrated data acquisition system consisting of an Instron Control 
Console, IBM PC/AT and software, a detail description and a list of procedures 
is in Appendix B. The object of the experiment was to test a current composite 
fiber type (used in many Navy aircraft) with a range of several gauge lengths 
and attempt to extract the underlying distribution of the fiber. The sample 
preparation is described in Appendix B. Each bundle sample was tested, the 
data stored, and quick analysis conducted to ensure no improper testing method 


were corrupting the data. 


A. BUNDLE TESTING 

The fiber material tested was made from a Hercules Magnamite high 
strength graphite, type AS-4 spool 145 in 3k bundles (i.e. there are 3,000 
filaments in a bundle). Several different gauge lengths were tested in order to 
analyze the effects of friction, slack, compliance on the fiber size. Size affect 
(Appendix A) should be able to correlate the location parameter beta between 
fiber lengths if weakest link physics is applicable. In all, nineteen bundle 
samples where tested. The test procedure was conducted as specified in 
Appendix B. The complete test log is listed in Table I. 

The bundle samples were handled with extreme care so as not to damage 
the any of the fibers. Moving the bundles from assembly bench to the storage 
bench and then to the testing machine presented the greatest hazard. The 


difficulty came when picking up the sample. This was done manually and 
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Sample # Gallp Oil/ Comments 


(cm) Dry 

090901 5.0 dry Good 

090902 0.5 dry No failure-slip 

090903 Os dry No failure-slip 

090904 Ons dry No failure-slip 

090905 0.5 dry No failure-slip 

090906 0.5 dry No failure-slip 

090907 5.0 dry Slip 

090908 5.0 dry Good 

090909 90 oil Good 

090910 Das dry Possible damage 

100901 aD dry Good 

100902 2.5 oi] Slip 

100903 2.5 oil Slip 

100904 50.0 dry Good - (friction) 

100905 50.0 oi] Good - (friction 
reduced) 

100906 020 oil Possible damage 

100907 25.0 01] Good 

100908 250 dry Good 





bending the bundles slightly was unavoidable. The long bundles (500 and 250 


mm) were less vulnerable to damage because of the relatively large aspect ratio 
but these were difficult to steady during transportation and tended to bounce 
around. 

Placing the sample in the grips also presented operational difficulties when 
tightening the upper grip with a socket wrench. An attempt was made to 
manually counter torque the grip, but this tended to be jerky at best. Also the 
grip pressure was difficult to control between samples. The alignment of the 
bundle in the grips was another concern. This was done by sight. If not done 
properly the bundle failure mechanism might be influenced. Manually reducing 
the experimental slack (not to be confused with bundle slack) was very touchy. 


When toggling the grips at the control console the load read out had to be 
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acutely monitored for the bundle loading must be kept to a minimum prior to 
testing. 

During testing a real time graphical display of the load vs displacement 
was used for monitoring the results, with this experimental irregularities in the 
several of the samples was detected. The bundles in the very short gauge 
lengths showed irregular load curves. Keying into this after several samples 
helped to identify a problem with the sample preparation. 

Friction was expected to affect the long gauge length (250 and 500 mm) 
and not affect the shorter gauge lengths (25 and 50 mm). To test the affect 
each gauge length tested a pair was treated with oil and a pair was left 
untreated (dry). Tests were run on samples of 500 mm, 250 mm, 50 mm, 


Eoemm and 5 mm. 


B. DATA REDUCTION 

After tests were complete the data files were decoded into ASCII and 
stored for processing. The data output from the Instron software is a three 
column file. The first column is the record number, the second is the 
displacement (inches) and the third is the load (Ibm). This is then processed 
through a program to convert the data to mm and kg and remove compliance. 
There are two data files output, one file with compliance removed and the other 
converted raw data (compliance not removed) [Ref. 5]. The output from the 
conversion program is then processed through an analysis program listed in 
Appendix C. This outputs the ECDF, the bundle Modulus E and the Slack 


region ECDF. The ECDF is then processed through a Weibull Maximum 


Likelihood Estimator program to define the shape parameter alpha and the 
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location parameter beta. ECDF is also processed through a linearizing program 
which outputs the Weibull probability plot data file. 
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V. RESULTS 


A. TEST OBSERVATIONS 

During testing several of the specimens were invalid due to experimental 
difficulties. Some of the samples had maximum strains noticeably beyond 
expected limits. One bundle was observed to have a much lower modulus than 
what was expected from equation C.5. Dissection of the samples with suspect 
loading curves revealed a lack of adhesive wetting resulting in a large 
percentage of the fibers slipping from the tab. The problem appeared most 
frequently in the small gauge length (GL) Samples. It is suspected that the 
adhesive which couples the bundle and tab was not wetting the internal fibers 
of the bundle, with a resulting loss of ability to carry the higher loads of the 
small GL. This phenomenon manifested itself two ways. If there was partial 
slipping after initiation of the failure region, a strain higher than the actual 
maximum strain was observed. If total slipping occurred in a large group of 
fibers the modulus decreased. After closer inspection of all the samples it was 
seen that many of the specimens had a small number fibers slipping in the tabs 
but the percentage was small enough not to effect the results. 

To avoid this problem the remaining tests were conducted in samples at or 
beyond 25 mm. The problem was less acute in the long GL because of the 
lower loads these long GL’s experience as a consequence of size effect. 

Two types of test were carried out, with oil and dry. This data was then 
reduced by categories of compliance removed and not removed. The compliance 
was estimated by testing samples of zero GL [Ref. 5]. The resulting curves 


were fitted and reduced to an analytical equation, displacement as a function of 
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load. This equation can then be used to subtract the compliance from the 
empirical data. The zero GL samples used the same coupling adhesive and 
possibly suffered from partial slipping. This causes the compliance curve 
coefficients to be over estimated and as a result displacements removed for 
compliance were too large. When the compliance was removed from the 
experimental data displacements in the lower range showed negative values 
indicating over estimations. 

Since the determination of compliance is slightly in error it is not used in 
the results discussion, though data reduction calculations were completed to 
indicate possible trends and affects. Table II shows the data reduction with 
removal, Table III is without compliance removal and Table IV, for comparison, 
is the result from single fiber tests [Ref. 6]. 

Table II. EXPERIMENTAL DATA RESULTS, COMPLIANCE REMOVED. 





Compliance Removed 


Test # Gage Modulus alpha Beta 
dry/oil length (kg/mm) 
(mm) (gm/e) (mm/mm) 

100901 25 109.76 915 3.0 0.0157 
dry 

090901 50 64.09 1068 Saal 0.0135 
dry 

100908 250 11.0 917 oo 0.00972 
dry 

100904 500 J436 905 3.035 0.0076 
dry 

090909 50 54.7 912 4.465 0.0145 
oil 

100907 250 11.18 932 3.964 0.0101 
oil 

100905 500 8 (es 959 3.66 0.0085 
oil 
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B. DATA COMPARISON 

Generally the results indicate only small deviations in the alpha parameter 
(Table III) for all GL. Between a oiled sample of the same GL only the 500 
mm GL had any noticeable change in the shape parameter, which had been 
expected. The location parameter was effected most by the oil treatment, beta’s 
increased by 10% in both the 50 and 500 mm GL. In 250 mm GL the alpha 
and the beta parameters seem to be relatively stable, meaning the oil treatment 
had little effect. 

It was expected that friction would affect the longer GL’s and have little 
or no affect on the shorter samples. The 500 mm test did show a considerable 
difference when oil was applied (Figure 5). 

The friction was expected to become influential only as the percentage of 
failed fibers increased, but the dry bundle shows a decreased load early, which 
can be seen in the ECDF plot (Figure 6). In fact the oil and dry sample data 
merge as the percentage of failures increases. The reduced friction made a 
difference in the resulting Weibull parameters, alpha increased from 3.085 to 
3.8. From single fiber testing alpha (Table IV) for AS-4 is approximately 4.11. 
Even with oil the bundle is affected, as seen by the sudden drops in the bundle 
load (Figure 5) or the sharp increase in the percentage of failures (Figure 6) at 
discrete points. 

Without friction the failure distributions are expected to be smooth. 
Decreasing the GL it is anticipated that the effects of friction diminishes. Tests 
at 250 mm the loading curve is seen to smooth out (Figure 7), providing a 
slightly more continuous plot. Friction had less effect on this GL though a 


slight decrease is noticed when the bundle is dry. 
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Table III. 


Compliance Not Removed 


Test # Gage Modulus Alpha Beta 
dry/oil length (kg/mm) (mm/mm) 
(mm) (gm/e) 
090901 30 45.6 761 ATS 0.01615 
dry 
100908 250 9.85 820 a07 0.0104 
dry 
100904 300 Dee 867 3.080 0.00781 
dry 
090909 50 AUS aS 4.89 0.0179 
oil 
100907 Zo0 a6 833 3.95 0.01041 
o1] 
100905 500 a0 917 3.8 0.00868 
oil 
SUMMARY (e0 Std if} Std. E 
OIL at Al 0.59 0.0162 COG 825 
DRY Selo 0.84 0.015! 0.00125 816 
Merged 4.09 0.66 OOTs7. 0:00 821 


DATA REDUCTION, COMPLIANCE NOT REMOVED. 


Beta 
sized to 
50 mm 

0.01615 

0.0154 


0.0137 


0.0143 
0.0154 


0.0152 


Std 


96 
Od 
10 





Table TV. RESULTS FROM SINGLE FIBER TEST AT 50 mm GL. 


Fiber Load Weibull Parameters 


Spool # Alpha Beta 
008 4.28 0.0175 
019 3.94 0.0172 
avernee — 4.11 —_ 0 0173 


Modulus = 900 gm/e 
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Figure 5. 500 mm gage length load curve comparison. 
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Figure 6. ECDF 500 mm comparison. 
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The ECDF plot (Figure 8) indicates friction had more influence than was 
revealed in the load curve. The upper tails again merge, indicating either the 
oil no longer decreases the friction or friction no longer affects the remaining 
fibers in this region of the ECDF. A slight change can be seen between the 
two samples in the lower and mid region but this could be variations between 
bundles and not whether the specimen was treated with oil, but it should not 
be dismissed because it is consistent with the other samples. 

A 50 mm GL test (Figure 9) showed results which had been anticipated in 
the longer lengths. Little change is observed as few fibers fail but the affects of 
friction become prominent as the percentage increases. 

The ECDF (Figure 10) indicates the effect of oil more clearly. The entire 
range of failure is shifted not only at the upper tail but the lower tail as well. 

Determining how well the underlying Weibull distribution is extracted from 
a bundle test and what influence the oil treatment has can be observed in 
Weibull plots. For the 500 mm sample (Figure 11 & 12) friction has definitely 
affected the distribution, but it is also observed that oil helped to reduce some 
of the deviation especially in the lower tail. 

As the gauge length is reduced it becomes apparent from the Weibull plots 
Mereures 13,14,15 & 16) for the 250 mm GL and the 50 mm GL that friction is 
not necessarily negligible, for the distribution is affected. This is clearly 
demonstrated in the 50 mm GL treated with oil (Figure 15) which shows the 


characteristic Weibull distribution distinctly as apposed to the dry sample. 
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Figure 7. 250 mm gage length load curve comparison. 
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Figure 9. 50 mm GL load comparison. 
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Figure 10. ECDF for 50 mm comparison. 


2d 





OS 18: 


(VL98/X)907 
GOO OS O- SO) ei ia 


1D Ww Oo0oS SO6O00!1 
LOld TINEAIM 
VIVG IWILNAWIdsdx4 


—————— 


OG t— 





Ore 


OS Ge 


00°- 00 s— tit 


OO 7— 


BOC 700. c— 


WOrC 


500 mm GL Weibull plot, oil treated. 


Figure 11. 
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Figure 12. 500 mm GL Weibull plot, dry. 
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Figure 13. 250 mm GL Wiebull plot, oil treated. 
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Figure 14. 250 mm GL Wiebull plot, dry. 
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Figure 15. 50 mm GL Wiebull plot, oil treated. 


32 


OS 'O 


(VLSE/X)901 
Ow oO OS @— Ol OM <a 


JO WW OG LO6060 
SO lap aes 
VIVQ IWLNAWIGddXx4J 


OS Lis 


CIC) a — 





OC 0am 


O0%- O0'%- 009- 008 
* 


000 


ODOC 





Figure 16. 50 mm GL Wiebull plot, dry. 
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VI. SUMMARY AND RECOMMENDATIONS 


A. SUMMARY 

The affect of friction is evidenced in all gauge lengths and is not always 
limited to the upper tail. Lubrication of the bundle reduces the effect of friction 
and should be used in all bundle failure testing. 

The compliance removal was not used in the data reduction due to the 
error in the compliance curve caused by adhesive coupling failure in the tabs. 
The compliance ended up being over estimated. This a side the trend of 
compliance removal on the data reduction (Table II) indicated a _ general 
reduction in the parameters (a,8). Estimating the compliance to be much less 
and extrapolating this thought, it can be assumed the value for the parameters 
presently determined would be slightly less. In a comparison of the single 
fibers parameters with the bundle parameters (Table JI), the bundle 
parameters were barely higher. Given the affect of compliance it can be 
concluded that once removed the bundle parameters will be indistinguishable 
from the single fiber testing and in fact given the quantity of the statistical 
base much more reliable. 

It was shown that friction effected all the gauge lengths but the longer 
gauge length (500 mm) expenenced the most disturbance. The 50 mm gauge 
lengths when treated with oil produced a very smooth distribution in the 
Weibull plot. The dry 50 mm sample produced a poor Weibull plot. The 50 
mm gauge length produce good results through out the entire range. The 500 
mm gauge length produced good results in the Weibull plot at the lower tail 


(Figure 7). Since the lower tail of the large gauge length are statistics of much 
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higher reliability than the 50 mm the linear statistics of the 500 mm can be 
pieced into the 50 mm statistical plot there by extending the range of the 
variable (strain). The upper tail of the 50 mm gauge length is not completely 
linear in the Weibull but this can be removed and statistics from a 5 mm test 
can replace it. This piece meal work on the Weibull plot can be done because 
of size effect. The slope alpha does not change but beta decreases in the 
Weibull as the gauge length increases. With this procedure friction and other 
contaminates can be removed and the underlying parameters accurately 
predicted. 

Through simulation a normal or uniform distribution of slack was 
determined to have very little effect on the distribution if the dO displacement 


is subtracted. 


B. RECOMMENDATIONS 
To ensure better results several areas of experimentation can be improved 
on. 
- The handling of the specimens should be reduced, a rigid carrying 
mechanism which can transport the bundle from assembly to storage 
to testing without any human handling or disturbing the fibers is 
necessary. 
- The manual mechanical grips should be replaced with a pneumatic device. 


- Qil treatment should be used on all samples. 


- Improve the tab coupling. 


The slack distribution requires more investigation than was addressed here. 
Simulation was a key element which helped to determine it’s effect on the 
distribution. An analysis algorithm is provided in this thesis which will help to 


determine the slack ECDF. To further study the question of slack several 
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bundle samples (large gauge length) dedicated solely for the determination of 
slack should be tested. A much larger data base is required to use the methods 
discussed in Appendix C. By decreasing the cross-head rate and increasing the 
acquisition rate a large data file can be attained. Once the distribution is fitted 
to a model, simulation can determine what affect it has on the bundle strength 


distribution. 
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VII. CONCLUSIONS 


The results achieved in this thesis conclude convincingly that the 
distribution of the fiber can be retrieved through bundle failure tests. The 
parameters extracted from the testing compared favorably with the single fiber 
testing. As a result, enormous characterization efficiency is realized; in fact 


with this technique, statistics to 10° are realized. 
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APPENDIX A 
SIZE EFFECT 


The determination of the scale parameter beta of a gage length, which is too 


short to preclude practical experimental implementation, can be accomplished by 


size effect of the weibull distribution function. 


Starting with the Weibull CDF for a fiber of defined gage length 1,: 


Pix) = Pf -exp Goa) 4) (A.1) 


R(x) = 1 - F(x) (A.2) 


For a second fiber of gage length 1, with, 


Lely Aone ily ent ta = nteecr a0 


and assuming the reliability R, is constant. 


R, = (R,\(RXR,)...(R,) = R, (A.3) 


substituting A.2 for R, 


R, = expl-(x/8,)")™ (A.4) 
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and the total reliability 


R, = exp{-(x/B,,)°") 


equating equation A.4 and A.5 


exp{-(x/B,)°"} = exp{-(x/B,)")™ 


taking the natural log A.6 reduces to 


(x/B)™ = m(x/B,)" 


again taking the natural log of A.7 


o.In(x/B,,) - o,ln{m(x/8,)} = 0 


pice C= & = a is not = 0 


BB. =m 
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(A.5) 


(A.6) 


(AC7) 


(A.8) 


(A.9) 


APPENDIX B 
EXPERIMENTAL PROCEDURES 


A. EXPERISNIENTAL SET UP 

The Test apparatus consisted of an INSTRON Model 4206 material tester, 
with a 50.0 kg load cell, connected to a 4200 Series Expanded Control Console. 
The Console was connected to a IBM PC/AT through com port 1 via a JEEE 
connection and converter. The data acquisition and material tester control was 
commanded by Instron Series [X Automated Material Testing System series 
4.01C software. 

The controlling software has a 5 page menu which requires setting prior to 
testing. Many of these settings are for Instron data reduction processing which 
are not used in this thesis but required values to operate. The settings crucial 
to testing and proper acquisition are cross-head speed, dala acquisition rate and 
gage length. Safety features include a maximum load and displacement setting 


to insure the load cell and testing model are not damaged. 


B. INSTRON CALIBRATION AND TESTING 
The Instron 4206 and control unit required a warm up time of one hour 
prior to operation to allow for stabilization. The load cell has two methods for 
calibration, an electronic calibration and a mechanical calibration. A mechanical 
calibration was performed for all tests. 
1. Mechanical Calibration 
Any inputs refer to the 4200 Sernes Expanded Control console. 


- Main power on/off/on - wait for diagnostic 


AO 


- Press Load Balance / Enter 

- Hand prescribe a 5.0 kg weight 

- Press Load Cal 

- Enter the weight / Enter 

- Remove the weight 

- Enter Load Balance / Enter 

- Re-hang the weight to verify correct calibration 
- Repeat if required 


2. Loading the Sample 


The procedures for sample placement into the test apparatus are as 


follows. 


- The Instron specimen grips are separated enough to 
allow the sample to hang from the top gnp without 
touching bottom grp. The tab of the specimen is 
firmly clasped with tweezers and gently lifted off 
resting bench allowing to hang free. 


- The samples top tab is placed in the top gmp, 
careful to align the bundle with the center, which 
has been measured and marked. After alignment the 
grip is tightened making sure not to twist or damage 
the specimen. 


- A load balance is performed at the control panel to 
remove the weight of the specimen. 


- The cross-head is then toggle down so as to place 
the bottom tab in the bottom grip then tighten. 
Note: after tightening the grips the specimen develops 
noticeable slack and the load cell measures a negative 
load. The specimen fibers are slightly compressed and 
need to be straightened. 


- The cross-head is toggle up until a load of 0.1 kg 
is indicated in the load readout window on the 
control panel. The load is then brought down to 
zero and the displacement reference zeroized. The 
procedure is repeated once to ensure the 
displacement returns to at or near zero. The 
specimen 1s ready for testing. 


3. Testing 
The test is initiated at the IBM PC/AT through the Instron software. 


After the specimen test parameters are set and loaded, prior to test initiation, 
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the IEEE port is enabled. The software then signals the system is ready to 
begin testing. 

As the cross-head displaces, load and displacement data are sent to the 
PC/AT at regular intervals of displacement, as prescribed by the initial software 
settings. 

A real time graphical output is available to monitor the testing. This 
allowed for a quick identification of bad tests. After the test was complete, the 
data file is converted from Instron System code to ASCII for latter data 


analysis. 


C. SAMPLE PREPARATION 
The samples were made from a Hercules Magnamite high strength graphite, 
type AS-4, spool 145. The bundle has 3000 fibers with a denier of .0057446 
gm/in.. 
1. Procedures 
A length of bundle is clamped at one end on a 3 meter aluminum track 
and weighted on the other. Half Copper tabs approximately 2.54 cm (width) x 
2.54 cm (length) are pre-positioned on the track at a specified gage length. 
Slots built into the track, center the bundle. An adhesive is applied at the tab, 
the bundle and the remaining half of the copper tab are securely clamped on to 
complete the tab. The bundle is severed between samples on the track and the 
adhesive is allowed to cure. 
Once dry the samples are manually removed from the track and placed 


on the sample bench. 


OO ———<« ee 


D. CALCULATION OF CROSS-HEAD SPEED 

The cross-head calculation is a function of the number of data points 
required, data acquisition rate, load cell and displacement gage tolerance. 

The time to bundle failure was selected to remain constant. This will 
provide a constant strain rate on the different gage lengths that were tested. 

The relationship between the cross-head, acquisition rate and recorded range 


is given in equation B.1. 
XH = (RRXAR)Y/DPTS (B.1) 


where: XH = cross-head speed 
RR = recorded range 
AR = acquisition rate 
DPTS = number of data points recorded 


The sampling interval is 
dena = Rhea OPIS = XH / AR (B.2) 


In order not to sample beyond the abilities of the acquisition system the 


following rule must be followed 
deltaD = Recording Tolerance (Bee) 


The Instron output tolerance for displacement was 0.00254 mm, therefore 
deltaD = 0.00254 mm. Using this as a guideline and the constant testing time 


chosen of 5 min.. Table IB lists the Instron test setting calculations. 
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Table IB... TEST SETTINGS. 


Gage Max Time Cross- Acq. 
length displ. head rate 

(mm) (mm) (min) (mm/s) (pts/s) 
500 25 5 8.33x107 6.6 
250 L275 5 4.17x10” 3.0 
50 Jigs 5 SSaxlo™ 3.3 
25 1.25 5 ATi txt Or on 
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APPENDIX C 


DATA ANALYSIS 


The analysis of the bundle failure load vs displacement curve can be divided 
into three distinct regions (see Figure 1C): 
- Slack region. 
- Linear region. 
- Non-linear region. 


Slack is created during specimen production and placement in the Material 


testing machine. As the bundle is measured, cut and tabs are placed on each 





LOAD CURVE 
REGIONS 


as) 


failure 








Figure 1C. Regions of the load curve. 


end of the sample some of the individual fibers become loose and effectively 
longer than other fibers. This produces an uneven loading across the fibers as 


the cross-head of the tensile testing machine strains the bundle. The i" fiber is 


not loaded until the displacement of the cross-head is equal to the respective 
slack of the i" fiber. 

When the cross-head displaces beyond the maximum slack, the loading curve 
becomes linear until the first failure. 

The amount of slack varies according to an undefined distribution. The 
slack distorts the failure region thereby disturbing the underlying sequence of 
the strength distribution. The effect slack has on the underlying shape 
parameter depends on the slack distribution and amount of slack relative to the 


total strain, this is demonstrated in Appendix D. 


A. SLACK DISTRIBUTION 


The slack in a bundle is dependent on the variation in the fiber lengths. 
L, = 1 + delL, (C.1) 


where: Lj = length of fiber j 
] = mean gage length 
delL, = slack of fiber j 


Since the delL, « | the length is approximated by 
L =] (C.2) 


J 


The load in the slack region can be defined by Hooke’s law. 


k-1 
P.(d,), = 2 E ((d, - dell,)(1 + delL,)} (C.3) 
jel 


k-1 
=~ZE d/l 2 k= 
j=l 


where: P.( ), = bundle load as each fiber k is loaded 
E = fiber modulus (gm/strain) 
d, = apparent displacement = d + delL 
d = fiber displacement 
n = number of fibers 
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PDF Plot 








Defining the loading slack in terms of the distnbution function (f,) of Ps at 


a displacement d, as shown in Figure 2C. P, can be write as 


Pd), 


(f(d,,)(d,5-d,,)n} E U(d-d,,){(d-d,,)/} (oy) 
where: d 


any displacement such that d > d,, 


U(d) = unit step function {=1 d 
{=0 d 


ams ee 
< 


da 
Summing the load (eqn. C.4) over the range of d, and defining the bundle 


modulus as a function of the fiber modulus 


E = nE/ (C.5) 


Pd) = E = [f,(d, dy j.. - dy J] Uld-d, ;) (C.6) 
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differentiating with respect to d for a continuous function 


(a 

6P/&(d) = E f-(d) did) = E Etc (Cae 
Jd, 

F(d) is the CDF of the slack. Assuming the empirical data approximates 


the continuous data. 


Posi. =: Weak (C.8) 


The derivative of the empirical data can be numerically solved with two 
methods, a numerical differentiation or a discretization. 
1, Numerical Differentiation 
Even interval Forward and Central difference methods are used due to 
the nature of the slack loading and the recording system. As the fibers load 
the slope is discontinuous as seen in Figure 3C. Several numerical methods 
exist which will produce’ varying 
accuracy. In trying to recover the 


underlying distribution of the slack Slack Resion 
| & 


(Pepresentation) 


the sampling rate and number of 
fibers must be considered (i.e., how 
many data points can be recorded 
relative to the number of fibers). 

If the sampling rate is fast 


and the number of total data points i 
« underlying loading 


. = : 5A: ' 
recorded in the slack region is 2-3 discrete data 


times that of the fiber, a 3 point Figure 3C. Slack region showing 


discrete data on the underlying loading 
slope. 
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forward difference method may yield accurate results due to the linearity 
between discrete data. If the relative sampling rate is slow the 3 Pt. forward 
difference may skew the distribution. To prevent this a central difference 
method is safe. 


3 pt. Central Difference method. 


lad = (-P ea 1 a hirelel (C.9.1) 


5 pt. Central Difference method. 


ee a (P 2° OP = Je! ot woven (Cn 052) 


ZL j4.y7OL 


where: hy eee 4 


2. Discretization 
The data in the slack region can be fit to a curve and a equation 
estimated. This equation can be differentiated analytically. The underlying 


data itself is not smooth so such a method is not recommended. 


B. LINEAR REGION 

The slope of the linear region is the bundle modulus. To calculate this 
slope a least squares method is very accurate, due in part to the natural 
linearity of the data. Before this method can be employed the region must be 
specifically defined. The point where the non-linear slack region ends is 


identified as dl and the end of the linear region is identified as d2 (the point of 
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the first failure). This is graphically displayed in Figure 4C. Between dl and 
d2 the least squares method calculates modulus EF. 

The points dl and d2 must be graphically identified and manually entered 
into the analysis program. To do this a plot of the load curve must be 


generated before the analysis can be conducted. 


LOAD CURVE 


ac) 





0 d, | 
Intercept point / qo od 


r4 
Figure 4C. The linear region. 


1. Least Squares Method 
Due to the noise, the discrete data 1 will deviate slightly from the 


underlying slope E. From Hooke’s law the expected relationship is 


T 
7] 
S 


P, = DV, (C.10) 


where: D, = d,-dl (C.10ne 


o0 





Figure 5C. Least Squares Method. 


But due to noise the load P deviates by an amount Pe, this is 


graphically represented in Figure 5. 


Pee— 7h - P. (Cit) 
Defining the least squares sum M 
M(E) = 2 {(ED, - P,)/20°}? (Cale) 


where: o = standard deviation 
Differentiating M with respect to E and setting it equal zero defines the 


minimum least squares. 


dM/dE = 0 


ol 


substituting in C.10.1 this reduces to 
E= YP(d,-di/2X(d,-d1) (C.13) 


2. Intercept of the Abscissas 
With the Modulus E the intercept of the abscissas can be calculated. 
q0 5 = dl =2(aie (C.14) 


C. FAILURE REGION 

This region is the key to retrieving the strength distribution. Using the 
slope of the bundle the failure statistics can be transferred to a ECDF space. 
From the ECDF the parameters of the distribution can be estimated. 

The ECDF is a plot of the percentage of failures versus a vanable. 
Typically this variable is the strength, this would be difficult to translate due to 
the nature of the bundle test. The displacement is constantly incremented, 
where as the load is a bundle load and can vary. Because of this the ECDF is 
plotted versus the displacement. ‘The displacement is proportional to the 
strengtli as defined by Hooke’s law and once the fiber modulus is identified the 
displacement can be translated to load. 

When a fiber fails, the bundle modulus E changes by one fiber strength. 
From Hooke’s law it is expected that the failures would occur along the 
modulus E;. E, being the modulus of the unfailed bundle which ideally is a 
multiple of the single fiber modulus proportional to the number of remaining 
fibers. Unfortunately this is not the case, as the bundle loads, friction and 
slack change the distribution causing the actual loading path to deviate from &, 


and follow along E,. 


oz 


To show this in more detail consider first the ideal underlying loading (no 
slack, friction or noise) of a bundle under constant strain. At a break in one 
fiber j, the load will drop until intercepting the slope of E,, shown in Figure 6C. 
This new slope emanates radially out from the intercept of the abscissas, (note: 
without slack the intercept should be the origin). If the data acquisition rate 
where high enough this failure could be clearly identified. Now add slack, this 
extends the _ failure 
displacement of some | 
fibers and the intercept || [| Failure 


of the abscissas by E 


) # fibers 
remaining 
— underlying 
loading 


after the initial failure 
to some point between 
the ongin and the dO ||, discrete 
point. The intercept data 

will shift a_ relatively 


small amount but will 





digress toward the 
Figure 6C. Failure loading path. 

origin. If friction is 
present the broken fibers will cling to the unbroken fibers and cause additional 
strain. Add noise to the acquisition system and the discrete data fluctuates. 
These factors combine to force the actual recorded data to deviate from the 
underlying bundle loading path. 

A bundle with n number of fibers will have n failures. If the data 


acquisition system can not record this many points in the failure displacement 


range than each statistic can not be recorded. With slack, friction and noise 
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even with a fast enough data acquisition system the failures might not be 
distinguishable. 

The percentage of failures in the bundle test can be defined by the ratio of 
the modulus of no failure to the modulus of failures. The reliability is defined 
by 

R, = E/E = P(d,)/E(d,,-d0)] (C.15) 


The percent failure is 1 minus reliability, the ECDF is defined by 


Fid,)=1-R (C.16) 


This relationship is graphically depicted in Figure 6C. The ECDF begins at 





Figure 7C. Transfer from loading curve to ECDF space. 
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d2, so the manual selection of this point must be done with care. 
Assuming all conditions ideal the ECDF would fall out as 
fee qa (Cla) 
But because conditions are not ideal this does not occur. What is important is 
the shape of the ECDF and whether it is close to the underlying CDF. The 
ability to distinguish each failure is unimportant. Another point to note is the 
lower tail needs to be clearly defined. Due to the relatively low failure rate in 
the lower tail and a constant recording rate this area of the curve has a high 
number of discrete data points per failure, so in essence the testing method 
takes care of itself. 
1. Upper Bound 

The data is recorded on even intervals of displacement not just at the 
failure points, this is graphically described in Figure 6. The discrete data most 
closely resembling the failure point is on the top edge of the loading curve. This 
is referred to as the Upper bound. When this is translated into the ECDF 
space the Upper bound is on the far 
edge of the plot as seen in Figure 8C. 
To retrieve the underlying CDF from 
the ECDF the Upper bound (Fu) 
must be extracted. 

Attempting to extract the 
Upper bound from a bundle of n 
fibers the most failure points that 
could be optimistically retrieved is n. 


With a large n this is not possible 





due to the acquisition rate, no Figure 8C. Close view of the ECDF. 
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matter, the shape and location is important and does not require a large 
number of points. 


Using the empirical ranking as a reference point 


F=k/m; k = 1,2,......m (C.18) 


where m is the number of points desired in the ECDF. F is used to define 
search bands of data where the upper and lower band is defined by 1/n. The 
value of F, which has the largest displacement (d,,) 1s assigned to Fu,. When 
no discrete data is in the search band F, Fu, 1s assigned the value of F at a 
displacement d, ,3. 

2. Maximum Likelihood Estimator. 

With the Upper bound extracted this data can be used to estimate the 
parameters of the distribution. The model selected is a Two Parameter Weibull 
Distnbution. The estimator chosen is the Maximum Likelihood Estimator 
(MLE) [Ref. #:p. 103] developed by R. A. Fisher. This method finds the 
maximum of the likelihood of the sample in the range of the vanable du and is 
a function of distribution parameters. For the Weibull model the method solves 
for the two parameters alpha and beta. The parameter beta is solved explicitly 
(Ref. 4] in terms of alpha. Alpha can not be solved explicitly and an iteration 


must be used. 


rare) 
ii 


[ Vm 2 due) * ~ “sivas (C13) 


{ © du.’ In(du,)AZ du,’) - © In(du,)} (C.20) 


R 
li 


o6 


An Atkinson method of iteration was chosen, initially alpha is assumed. From 


this successive alpha’s are calculated using equation C.20. 


> 


o, = a, - (a, - a,)*o, + 20, + o) (Cz 


where: «a = initial start value 


oO 


Oh, = f,( a, ) 
a, = £,(0) 


This equation is iterated for a’, until o’, = q,. 
1. Linearized ECDF 
The data is linearized for a Weibull probability plot to visually check 
the conformity to this model. To linearize the ECDF start with the Weibull 


equation 
F = 1 - exp(-P/8)* 
subtracting 1, multiplying by -1 and taking the natural log. 
In(1-F) = (-P/B)* 
again taking the natural log 
In[-InQ1-F)] = odn(P/f) (e227) 
C.22 is a linear line with a slope of alpha. 


D. DATA OUTPUT 
As the analysis program calculates certain data file are output, these files 
are: 
- EXPER.OUT, the converted experimental output displacement, strain 
and load. 


- SLACK.OUT, the slack region. 
- EXMS.OUT, the load curve minus slack region. 


O7 


- UBFECDF.OUT, the Upper bound ECDF. 
- FECDF.OUT, the ECDF. 
= SECDFOUT, slack BODE 
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APPENDIX D 


A. SIMULATION 

The simulation of a bundle test has two purposes. First by analyzing in 
depth the mechanisms which govern the failure of the fibers in a bundle, a 
better understanding of how to extract the characteristic properties is attained. 
The properties desired are the distribution parameters. 

Secondly a simulation gives us the ability to know what the underlying 
parameters are and thereby proof the analysis data reduction algorithm. Also 
the variables which influence the shape of the failure distribution can be varied 
in simulation and the degree of induced affect observed. 

The bundle loading curve can be divided into three regions; 

- Slack 
- Linear 
- Failure 

Initially the simulation requires input of the vanable estimates, this is 
done through an input file. 

1. Initial Inputs 

The underlying parameters define the strength distribution for a gage 
length. The distribution chosen is the Weibull model. The parameters are 


alpha and beta, the number of fibers per bundle is n. The fiber modulus (E) is 


in gm/mm’, from Hooke’s law 


o=Ee (D.1) 
where: oO = stress 
€ = strain 


og 


Substituting the definition of stress and strain for fiber i 


P/A = E 8/L. (D.2) 


where: P = load 
A = area = constant 
L = length of fiber 
5 = change in length L 


rearranging D.2 and letting E, = EA (gm/e) 


P= Reo (D.3) 


2. Slack and Strength Distributions 
The slack is caused by the variation in fiber length L. If L were 
constant no slack would be evident. The slack is expected to have a definable 
distribution, but this is unknown and could vary from bundle to bundle 
therefore several models are simulated. Two models selected where 


- Two parameter Weibull 
- Uniform. 


Although physical consideration favors the slack model to be normal, 
the Weibull model was chosen for it’s analytical convenience. The Weibull is 
very flexible and can characterize a wide range of shapes. For an alpha of 3.5 
the distribution is approximately normal, other than 3.5 the shape of 
distribution is skewed negative or positive. Beta identifies the central tendency 
of the distribution, its relatable to the mean. From equation 2.1 the Weibull 
CDF is defined as 


F(x;0,8) = 1 - exp{-(x/R)* 
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rearranging in terms of the variable x, which can represent slack (delL,) or 


strength P.. 


x, = exp{ [In(-In(1-F,) + odn(B)] / a) (D.4) 
O< FL, <1 


The uniform model is a simple two variable distribution with a lower 
(u,) and upper (u,) range. Since the slack is initiated at zero the lower range is 


zero. 


F (delL;u,,u.) = (delb - u,)Au, - u,) (D.5) 
Oe ae | 
solving for delL, from D.5 and setting u, = 0 


delL, = u, F (D.6) 


3. Data Generation and Continuity 
A uniform random number generator creates an array of numbers 

between O and 1 and stores it until recalled. 
Recalling equation D.3 and substituting in C.1 and solving for the 


displacement (8). 
ona + dell.) /-E, (D.7) 


The bundle is put under a constant strain rate, so the discrete 


variable will be displacement. As the bundle is strained the filaments will fail 
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in order of the underlying strength or strain (neglecting slack and friction). For 
a constant filament length the strain is proportional to the displacement and 
the fiber will fail at some underlying displacement. If slack is present in the 
bundle each fiber will not fail at the underlying displacement because the 
displacement is no longer proportional to strain but a function of both slack and 
stain. 

The displacement observed by the Instron material testing machine is 


not 6. The machine records an apparent displacement (d,) 


d,, = 5, + delL, (D.8) 


This will affect the fiber failure order, every other variable is then 
keyed to the order (n,) of d, ,. 
4. Slack Loading 
As the bundle is strained, the fibers are loaded by a progressive 
summation of each fiber in order of their slack distribution. The loading 
function in the slack region can be defined as a function of the displacement (6,) 


each fiber experiences after slack is released. 


je 
Peay 2 Een 4 =laZee ae (is 
i=l 
j-l 
P,,= YE, ((delL,, - dell, (1 + delLy;)] (D.10) 
j=] 
where: delL,, = order slack 
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This relationship can be represented graphically in Figure 1D. The last load 
calculated from equation D.10 will occur at delL, ,. 





Figure 1D. Slack region load simulation. 


4. Failure loading 
The order of each fiber failure is determined by the ordered strength 
distribution (assuming no slack). With slack the order of failure will deviate 
slightly. The failure load is a function of the modulus of the bundle and the 
displacement of each fiber. The bundle failure load (P,) can be defined from 


equation D.2 as 


yy he 10, (D.11) 


i=) 


substituting in C.1 and D.8 


P,,= LE [(d,, - dell, + delL] (D.12) 


1S 
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5. Discrete Data 
The material testing machine measures and transmits data on even 
intervals of displacement. To simulate this the defined loading points 
(P,,delL,,) and (P,,d,) will be set limits for a linear interpolation. The data is 
interpolated on intervals of deld, which is set by the cross head speed. 


Defining the discrete apparent displacement 


Ce dn dca (Dit 


The interpolation algonthm to define the discrete data from 
continuous is the same for both slack and linear region. 
The discrete load P in the slack region is a function of P., the ordered 


slack (delL,) and the apparent discrete displacement 


P, = ((P,,,, - P, MdelLy,., - delLy Id, , (D.14) 


a )el 


The discrete data is defined this way until the first underlying 
displacement d, .. 
The slope of the linear region intercepts the abscissas at a point dO. 


The slope of the linear region is the bundle modulus E (see Figure 2C). 


E = (P,, - P, Ad, , - dell.) (Dita) 


The slope can be solve using E 


d0=dlL-P JE (D.16) 
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LOAD CURVE 





Figure 2C. Linear calculations. 

Subsequent discrete data points can be created by a regressive 
process. At each fiber loss a new slope emanates radial outward from dO. 
Between breaks, to define the discrete data, the slope is defined by the failure 


load P,, the displacement d, and dO 


P, = [P, Ad, ; - dO)] (d,; -d0) (D.17) 


where: Cede —le2.....,T1 


6. Noise 
Noise was simulated simply by setting a tolerance level and using the 
random number generator to produce a set of random numbers between -0.5 


and 0.5. This data was added to each discrete load P. 
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The input file to the simulation has units of kg,mm and the 
expernmental output from the Instron Software is in lbm,in. A conversion of the 


data is required, this is accomplished through a separate program 


B. SIMULATION INTERPRETATION 
The simulation was written in fortran code, listed in Appendix F. Multiple 
simulations were run varying the parameters slack and noise to see the effect. 
Simulation data was analyzed via the algorithm as represented in 


Appendix C. The output is the bundle modulus E 


=—) 


the shape parameter of 
distribution Alpha, and the location parameter Beta (mm/mm). 

The input parameters were alpha = 4.0, beta = 0.16 kg and a 50.0 mm 
Gage length (GL). The noise was varied from zero to about 0.1% fluctuation of 
the maximum load. The maximum slack input was 3% of the gage length and 
a normal slack distribution. The seed used for the random number generator 
was held constant so the underlying strength distribution was constant. The 
seed for the slack region and noise was not held constant. The results of 
multiple simulation test is in Table IC. A value of 4.134 for the shape 
parameter was retrieved consistently (same seed), as expect no deviation for the 
same data. The increase in slack drove the parameter down. The noise did not 
have any effect on the parameters. Figure 3C show the loading curve with 
slack (minus dO) and no slack. 

No noticeable change in the parameter output nor in the ECDF plot 
(Figure 4C). 
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Table IC. SIMULATION 


PROOF TESTING RESULTS. 


Maximum Slack alpha 


o% GL 
(8 = 0/2) 


0 
5) 
10 


4.134 
4.010 


4.099 
4.095 
4.091 
4.089 
4.088 
4.075 
4.076 


OBSERVATIONS AND 


beta 
(kg) 
0.016 
0.016 
0.016 


The numerical derivative of the slack region is very sensitive to the any 


noise. 


This can be seen in the plot of the slack ECDF Figure 5C. 


In order to 


retrieve the slack distribution the recording system must have a low tolerance 


load cell. 
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Figure 3C. 
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Figure 4C. The simulation ECDF for slack and no slack. 
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Figure 5C. Simulation of the slack ECDF. Weibull approximation of Normal. 
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APPENDIX E 


I. DATA ANALYSIS CODE 
A. ANALYSIS PROGRAM 
1. Data Input 
The program requires a formatted data file called EXPER.OUT. The 

data is read in and prompts the user for several variables. These are listed 
below. 

- What is the gage length? 

- How many fibers in the bundle? 


- Which Numerical differentiation is desired? 
3,5,7 pt. forward diff. or 3,5 pt. central diff. 


- The end of the slack region and the beginning of the 
failure region in displacement (mm) (this requires a graphical 
estimate). 


- How many points in the input data file? 
- The acquisition tolerance of the displacement? (this is for 


any filtering of data points in the event the sampling rate 
is too high and repeated data is evidenced. 


- The modulus is calculated and displayed and then ask if this 
value needs to be change. This is in the event the modulus 
calculation in correct a manual entry can be made. (Be Careful 
not to corrupt the true data). 


2. Data Output 
The data out put is to the screen and to a file. The following is a 
list of the file outputs. 
- EXPER.OUT, the input file is output with an added column of 


strain 
- SLACK.OUT, the slack region (load vs displacement) 
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- EXMS.OUT, the load curve minus slack region 
- UBFECDF.OUT, the Upper bound ECDF 
- FECDF.OUT, the ECDF 
- SECDF.OUT, the slack region ECDF 
B. MAXIMUM LIKELIHOOD ESTIMATOR 
The MLE program reads one of the ECDF data files, users choice, the first 
line in the input file is a utility line defining the file size, the number of 
filaments and the bundle modulus E, the remaining lines are data. 
1. Data Input 
The input other than the read file is prompted from the screen. The 
following is a list of requested information. 
- How many points for the Expected Ranking (The MLE requires a 
probability ranking) 
First approximation for alpha 
2. Data Output 


The parameter estimations alpha and beta (mm/mm) are output to the 


screen. 


C. LINEARIZED PLOT 
This program takes the ECDF and does a Weibull linearizing (equation 
C20). 
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Il. FORTRAN CODE 
A. ANALYSIS PROGRAM 


PROGRAM ANALY 
BUNDLE FAILURE ANALYSIS 


Outline of program: 
Read in the data points 
NIA - Number In Array 
2K KK KK 2K IK KK KK KK 2 KK KK 2K ok 2 9 KKK 2 9 ok KK 2 Ko ok oi 2 ok 2 2 OK oi Ko OK 2 ok KK KK KK KK KKK 
4 Find the 
-linear region 
. -slope of the linear region = EB (E BAR) 
- -intercept point of the linear region with 
* 
os 
* 


%* & &€ He KH OF 


the abscissas = DO (d zero) 
-Maximum slack or minimum linear displacement - DP1 


-Maximum displacement of the linear region - DP2 
AK 2 KK 9K KK KK Koo 2K KK 2K KK 2K KK KK 3 KK KK KK KK KK KK KK Ko 2 KK KE Koo 2 2 KK KR 2 KK OK oR kK ok 


* 


* Input the range of the linear region by manually retrieving 
* the defining points DP1 and DPz2. 


* 
2 KK 2 2 2K 2 2K 2K 2k 2 2 KKK 2K 2K oi 2 OK 2 2 ok 2 ok 2 2 2 2 2K KK ok 2 2 KK 2 ok ok 2K 2K i 2 OK 2K 3 3 2 KK ok 2k OK ok 


* OUTPUT: SLACK.OUT  - slack region. 


‘a SECDF.OUT - slack ECDF. 

. EXMS.OUT _ - experimental data minus slack region. 

ss FECDF.OUT - failure region ECDF. 

ii UBFECDF.OUT - upper bound of the failure region ECDF. 
* 


INTEGER NIA,NISA 

PARAMETER(NIA=6000,NISA=1000) 

REAL D(NIA),P(NIA),PP(NIA),H,LP(NIA), LD(NIA),FS(NIA),F(NIA) 
+,FUB(NIA), DF(NIA),DFUB(NIA),PF(NIA),PFUB(NIA), TEMP(NIA),GL 
INTEGER COUNT,NIDP1,NIDP2,NI(NIA),SLCT 

CHARACTER*1, Q1,Q2 


* 
2 2 2 2 2 2 2 2 2 2 2K 2 2 2 Ko Ko KK OK 2 2 KK 2 Ko kK 2 oR 2K Koo 2 KE i 2 2k 2 kK KK KKK ok Ko ok 2 ok ok 


* Input analysis control parameters. 
* 


36 2K 2K 2 2 2 2 OK 2K OK KR oo 2 2 KE 2K KK Ko KK oo Ko KK Ko Ko Koo oo oo oo oo 2 2 ok ok 2K ok 


PRINT *, "What is the gage length? (mm)’ 
READ *, GL 


PRINT *, ’ How many filements” 

READ *, N 

PRINT *, ’ Which method do you require for reduction of 
PRINT *,’ the slack region to the ecdf space?’ 
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PRINT * forward diff. method.’ 


ae 1) 3325 
PRINT *, ’ 2) 5 PT forward diff. method.’ 
PRINT =~ 3) 7 PT forward diff. method.’ 
PRINT *,” 4) 3 PT central diff. method.’ 
PRINT 5) 5 PT central diff. method.’ 
READ *, SLCT 


* 
A Re 2 a 2 2 ae oe a 2 2 ok ok 2 kc 2 ak ok 2 a 2 a ok a 2 ok ok 2k 2K ok oe a 2 2 ok ok 2 2c 2 2 a ok ie a ok 2k ak ae ok ok 2k 2k 2k ok ok 2 2 2k a Kok 
* Read the range of the linear displacement. 

* 

BE 2 ie 2 ate 2 ae ok a ake 2 2 2k ae oe a 2 2 2 ae a 2 2 a 2k a ok 2k ak a a ie aK ae ok 2c 2 ak a ok 2k ok ok 2k oe ok ok a 2k 2k a 2k 2 2 2 2 a 2 2k kK ok 2K ok 
* 


PRINT *, "INPUT DP1 (the maximum slack)’ 

READ *, DP1 

PRINT *, "INPUT DP2 (the maximum linear displacement)’ 
READ *, DP2 


~-S PEPE LE LES SSS FTES TEES TTS FF SST TSS TTT ST TST ST TTT TS TT TS TF TT TT TT TS tS St SS 


* Find the how many points are in the data file. 
* 


AE 2 2 2 Ee 2 2K 2 2 oie 2c 2 2c of 2c aE a ofc 2 af 2 a OE i 2 2K of 2K 2 oo oie 2 a 2 a oe 2 2 a 2 2 a 2 a a 2 oo ok 2 2 oi ok a OK ok 2 OK OK ok 


PRINT*, "How many points in the data file?’ 
READ *, COUNT 
22S SSE SST ESSE SLE E SESE LSS ESE SESS SSE STE SE SELES SEES ESSERE SESS TS TSS ST FSS SS SF 
* Find out what the tolerance of the machine is. 
* (displacement) 
(eter et SELES SS FTL ESTES LESTE SEES E LE TSE LEE SESESEEESIESLSEEELESLSEEESEIEST TEESE LESS +25 + F 3 
PRINT *, What is the tolerance of the displacement output’ 
PRINT *, data? This is the limit between data displacement’ 
PRINT *, Input in mm.’ 
READ *, TOLD 
PRINT *, "How many points in the empirical rank?’ 
READ *, NER 
144 PRINT *, "Do you want the load in Kg (K) or Newtons (N)?’ 
REAIX(*,303) Q2 
IF(Q2 .EQ. ’K’ .OR. Q2 .EQ. ’N’) THEN 
PRIN S27? 
ELSE 
PRINT *, Try again in UPPER CASE K OR N’ 
GOTO 144 
ENDIF 


* 
PS tS S SITS SS TLL SL SL SESS SLES SL SSS SSL SS SSS SSS SSS SSS SS TS TST SETS FS TTS T STS TS TSF FS 


* Read in the experimental results. 
* 


PSP SSF LS FSF S SSS FS FSF SS LES SS SST LSS STS SS SST SF FSS TTS FT FFF TTT SS SS SSS SS SS tS 
* 
OPEN(UNIT=11,FILE=7EXPER.OUT,STATUS=’OLD’) 
x 
DO 10 IJ=1,COUNT 
REAIX(11,110) N11), TEMP(1),P(1) 
10 CONTINUE 
REWINT(11) 
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id 


CONVER = 1.0 

IF(Q2 .EQ. ’N’) THEN 

CONVER = 9.8062 

ENDIF 

DO 44 I=1,COUNT 

P(I) = P(I)*CONVER 

WRITE(11,1001) NI(1), TEMP(),P(1), TEMP(1I)/GL 
CONTINUE 


FO OOOO domi aR ok a aC i a aca ak akc ak ak ak ake ak a akc ake a ak afc afc a 


* Filter out any points with difference in displacement interval 
* less than the machine tolerance TOLD. 
* 


OOO OR a Rm ok ogg icici ak ak ak ak akc ak ak ak ak akc ak 


30 


* 


FIND COUNT (NI) OF DP1 AND DP2 
* x * x * 


J=1 

D(1) = TEMP(1) 

DO 11 I=2,COUNT 

DIFF = TEMP(I)-TEMP(I-1) 
IF(DIFF .LT. TOLD) THEN 
PRINT *, “SKIPPED POINT ’,I 


ELSE 
J=J+1 
D(J) = TEMP(I) 
Pi) = Pi) 
ENDIF 
CONTINUE 
COUNT = J 
* * * * * * * 


I=1 
IF(D(I) .LT. DP1) THEN 
I=I+1 
GOTO 20 
ELSE 
NIDP1 = I 
DP1 = IKI) 
ENDIF 
IF(D(I) .LT. DP2) THEN 
I=I+1 
GOTO 30 
ELSE 
NIDP2 = I-1 
DP2 = IXI-1) 
ENDIF 


We Re ee ie oe oe 2K ie te 2 2 2 ie 2 a i 2K i 2 a i 2 a ok 2K oc ic a a oe a oe i 2 ko 2 ko 2k 2k 2c ok 2 ok 2K 2 oc 2k oe ok ke ok ok 


* OUTPUT THE SLACK REGION. 
* 


ARK OK oi a ik ak ak kk ak ak ak ke ok ke ac ak 2c oc ok ok ie ote oc a 2 ak oc 


* 
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OPEN(UNIT=13,FILE=SLACK.OUT,STATUS=UNKNOWN’) 
DO 500 I = 1,NIDP1 
WRITE(13,110) I,D(1),P(1) 
500 CONTINUE 
CLOSE(13) 


* 
AK A KK KK 2K 2K 2K 2K 2K KK 2K 0K KKK OK KK KKK KKK OK KK KK KK KKK KKK KK KK KK KK KK KK KKK KK KKK 


* LEAST SQUARES TO DETERMINE THE SLOPE BETWEEN DP1 AND 
DP2. 


* 
AE A a AR 2 2K Ko ok a ok a a oe 2k 2k 2 a a a a a a a OK 2 2 KK KK KK KK 2K KK ok 2K 2K 2k ok 2K Ko KK ok ok ok ok ok KK KKK KK KK 
£ 


0.0 
0.0 
I=(NIDP1+1),NIDP2 
LP(I) = P(I) - PONIDP1) 
LIX]) = DD) - DP1 
TOP = TOP + LP(I)*LD(1) 
BOT = BOT + LIXI)**2 
40 CONTINUE 
EB = TOP/BOT 
+ * 


TOP 
BOT 
DO 4 


On il 


x 


* Find the intercept point do 
x* * a = 


+ ak ae ye ak 3k 


x 


DO = DP1 - P(NIDP1)/EB 
PRINT *, "Do = ’,DO 
PRING *,”° 
PRINT *, The Modulus calculated from the load curve’ 
PRINT * 1s -EB 
PRINT *,’’ 
PRINT *, "Would you like to change this (Y/N)? 
READ(*,303) Q1 
IF(Q1 .EQ. 'Y .OR. Q1 .EQ. ’y’) THEN 
PRINT *, Input the new modulus’ 
READ *, EB 
ENDIF 
303 FORMAT\(A1) 
* 


*x 
2K aK 3K 2K 2K 2 2 oe 2K 2K ok 2K a 2K ok 2K ok 2K 2K 2K 2K 2K 2k KK ok 2K 2K 2 ok ok oo oo oo ook ok ok KK KK KK KKK KKK KKK KKK KKKKKE 


* OUTPUT THE EXPERIMENTAL DATA MINUS SLACK. 
x 


ae ae a ak 2 a ae a a 2c 2 a a 2 2 2 a 2 a 2 ok a ok 2 ok 2k 2k a 2k a 2k ok a 2k ok a 2k a ke 2k a ak ak 2k ok ak ok 2 2 a a 2k 2 oR ok ok ok ok OK KK KK KK KK 


OPEN(UNIT=13,FILE=EXMS.OUT ,STATUS=UNKNOWN’) 
J=1 
WRITE(13,110) J,0.0,0.0 
DO 550 I = NIDP1,COUNT 
J=J+1 
WRITE(13,110) J,(DC1)-D0),P(1) 
550 CONTINUE 
* 
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2 2 2 2 2k ok a ok i aK 2 2 ok ok 2 2K ok 2 oe 2K ok oe a oie ok 2 2 a 2c 2 aK oe ok 2 EK oe 2 2 oC 2 2 2 2 2 2 2 2 2 2K 2 2 a 2K ok ok ok ok ok 


* NUMERICAL DIFFERENTIATION OF THE SLACK REGION 
i (using a forward and central difference methods). 
3 


2 9 2 2 2 2 2 2 2 2 ok ok 2k 2k 2 2K 9k 2 ok ok ok ok 2 2 ok 2K ok 2 2 a 2 2 ok ok ok ok ok 2 2 ok ok 2 2 oi oo of 2 2 a ok a ok 2 2 2 oo Ko ok KK ok ok ok 
* 


IF(SLCT .EQ. 1) GOTO 103 

IF(SLCT .EQ. 2) GOTO 105 

IF(SLCT .EQ. 3) GOTO 107 

IF(SLCT .EQ. 4) GOTO 130 

IF(SLCT .EQ. 5) GOTO 150 
* 


2 2K 2 2K 2 2 2 2 2K ooo 2K 2 2 2K ok 2K 2K 2K 2K ok ok 2k ok OK 2k 2 KK KK 
* 3 pt forward diff. method. 

2K 2 2 2k 2 2 2K 2k 2 2 ok 2K 2 ok 2 2 2 2 ok ok 2K 2K a ok 2k 2k ok ok ok 2K ok KK 
aS 


103 DO 3 I=1,NIDP1 
* if out of the slack region leave the loop 
IF (D(I) .GT. DP1) GOTO 51 

H = ABS(D(I+1)-D{1)) 
* 3 point forward difference numerical differentiation. 
* Generate Fs(da) CDF 

PP(I) = (-3*P(I)+4*P(1+1)-P(1+2))(2*H) 

FS(I) = PP(I)/EB 
3 CONTINUE 

GOTO 51 


AK 2 3K 2K 2K 2K ok KK 2k 2K KOK kK ok aK EK 2 OK 2k Ko oi Ok ok ok OK ok ok ok ok ok 


* 5 pt forward difference method. 
2K AK 2K 2K 2K aK 2 2K 2K 2 ok 2K KK 2K ok 2K ok 2K 2K 2K OK ok OK 2K OK 2K 2K OK OK OK KK 
105 DO 5 I=1,NIDP1 
* if out of the slack region leave the loop 
IF (DUI) .GT. DP1) GOTO 51 
H = ABS(D(I+1)-DXI)) 
* 5 point forward difference numerical differentiation. 
* Generate Fs(da) CDF 
PP(I) = (-25*P(I)+48*P(I+1)-36*P(1+2)+16*P(1+3)- 
C3*P(I+4))X12*H) 
FS(I) = PP(D/EB 
5 CONTINUE 


AK A A A 2 I KK 2 kK EK OK 2k KK oe 2 2 oo 2 i oi oi ok ok 2K ok ok 


* 7 pt forward difference method 
2K 2 2 2K KK 2K 2 2 KK 2K 2 KK ok 2K 2K OK 2K 2K 2K 2K 2K OK ok OK KK OK OK OK OK 
107 DO 7 I=1,NIDP1 
* if out of the slack region leave the loop 
IF (D(I) .GT. DP1) GOTO 51 
H = ABS(D(I+1)-D(D)) 
* 7 point forward difference numerical differentiation. 
* Generate Fs(da) CDF 
PP(I) = (-147*P(1)+360*P(I+1)-450*P(1+2)+ 
+400*P(1+3)-225*P(1+4)+72*P(1+5)-10*P(I+6))(60*H) 
FSI) = PP(ID/EB 
7 CONTINUE 
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2 KK KK i KK KKK KKK KKK KKK KK KKK KK KKK KK 
* 3 pt central difference method 

AE 2K oe oe 2k 2 ok a 2k a oe 2 a a oo aK ok a KK KK ok OK OK KK 
* 


130 I=1 
H = ABS(DX1I+1)-D(])) 
PP(I) = (-3*P(1)+4*P(I+1}P(1I+3))/(2*H) 
FS(I) = PP(I)/EB 
DO 33 I=2,NIDP1 
* if out of the slack region leave the loop 
IF (DI) .GT. DP1) GOTO 51 
H = ABS(D(I+1}D())) 
PP(I) = (-P(I-1)}+P(I+1)/(2*H) 
* Generate Fs(da) CDF 
FS(I) = PP(D/EB 
33 CONTINUE 
JK 2 KK i oo KK i KK Ko KK KK KK KKK KOK KK 
* 5 pt central difference method 
3K OK 2 2 ie i oe i 2K oe kK ok oe ok ok ook ok ok ok kok KKK KK KK KKK 
* 


150 = 
* Use a 3 pt forward difference method for first couple pts. 
H = ABS(DX(1+1)-D(1)) 
PP(I) = (-3*P(1)+4*P(I+1)-P(1+3))(2*H) 
FS(I) = PP(YEB 
l=2 
H = ABS(DX(1I+1)-D(1)) 
PP(I) = (-3*P(1)+4*P(I+1)-P(1+3))2*H) 
FS(I) = PP(I)/EB 
DO 55 I=3,NIDP1 
* if out of the slack region leave the loop 
IF (D(I) .GT. DP1) GOTO 51 
H = ABS(DX(I+1)-D(1)) 
PP(I) = (P(I-2)-8*P(I-1)+8*P(1+1)-P(I1+2))(12*H) 
* Generate Fs(da) CDF 
FS(I) = PP(D/EB 
55 CONTINUE 


FE 2 2 2 KK KK OK KK RK KK KK KEK KKK KKK 


51 CONTINUE 


AE A A A OR OK EE 2 2 2 2 2 2 oo eo ei ee 2 oe oi 2 KK KK KK KK 


* FAILURE REGION 
* 


"eR KR ko Rk kk ok Kk KK kk kok kkk a ko oo ok ko ok ke ok kok ok ke ok 2c 2k ie 2 ok ko oo io oo oo KK ok KK 


J=0 
DO 60 I=NIDP2,COUNT 
J=J+1 
F(J) = 1 - PUMEB*(DXI)-D0)) 
DF(J) = D(I)-DO 
60 CONTINUE 
= J 


Pet SS SF FS SS SST TS SS STFS FS SSS SSS SF SS PSS SS FSS FS FSF FLT FPS S FFF TS SS SS SS St tS 


* EXTRACT THE UPPER BOUND. 
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* 


KK 3 2 2k 2k 2 2k 2K Ko 2k 2 i i 2k 2 ie 2k 2c 2c 2c 2k 2a 2 2k oi i 2k 2k 2c a 2K 2k 2k a 2c 2k ik 2c a 2k ik ic ic ic 2c ik 3K 2h 2c 2h ik 2k 2k 2k 2k 2 aK 2k 2k 2k 2c 2k 2c 2k 2k 


90 
91 


% 


93 


94 


92 


DEL = 1.0/NER 
K=NUMF 
I=NER 
DFUB(I+1) = DF(K) 
K=K-1 
I=I-1 
FL = 1.0*I/NER - DEL/2 
FU = 1.0*I/NER + DEL/2 
PRINT *,7FL F(K) FU’ 
PRINT *, FL,’ < ’,F(K),’ < ’,FU 
IF((K .LT. 1) .OR. (I .LT. 0)) GOTO 92 
IF((F(K) .LT. FU) .AND. (F(K) .GE. FL)) THEN 
DFUB(I+1) = DF(K) 
J=I-1 
FL = 1.0*I/NER - DEL/2 
FU = 1.0*I//NER + DEL/2 


K=K-1 
IF(K .LT. 1) GOTO 92 
GOTO 91 
ELSEIF(F(K) .LT. FL) THEN 
ae ok 3 a * aS so 
check to see how many it skips 
x + x x x ME 
Jl =1 
IF(F(K) .LT. (1.0*(1-J1/NER - DEL/2)) THEN 
Jl=J1+1 
GOTO 93 
ENDIF 


DO 94 J1sI,(I-J1+1),-1 
IF(11+1 .LT. 1) GOTO 92 
FUB(I1+1) = 1.0*I1/NER 
DFUB(TI1+1) = DF(K) 
CONTINUE 
I=I1 
FUB(I+1)=F(K) 
DFUB(I+1)=DF(K) 
I=I-1 
FL = 1.0*I/NER - DEL/2 
FU = 1.0*I/NER + DEL/2 
K=K-1 
GOTO 90 
ELSE 
K=K-1 
IF(K .LT. 1) GOTO 92 
GOTO 90 
ENDIF 
CONTINUE 


PRINT *, "OUT OF LOOP WITH K AND I’ 
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PRINT *, Keene la 
NUMEF = [+2 


2 ae 2 2K 2 2 2K 2 fc oe ke 2 kc 2 akc oe oc ae ke a 2 ake 2 2 2K a 2K oC 2K 2K oe 2 a 2K 2 2K 2 OK 2K KO OK a OK OK EK 


* WRITE THE DATA 
* 


OR AK 2K 2K 2 2 oe 2K 2K 2 2 aK oe oe 2K oe oe oe 2K 2K 2K 2K 2K 2K 2K 2K 2 2 eo oe ee eC oe oe 2K 2K OK 2K oe eK OK KE OK KK 


290 


70 


80 


OPEN (UNIT=11,FILE=UBFECDF.OUT,STATUS=UNKNOWN’) 
MSUM = NER+1-NUMEF 
WRITE(11,1010) MSUM,N,EB 


DO 250 I=NUMEF,NER 
WRITE(11,110) I,DFUB(I),FUB(T) 
CONTINUE 
CLOSE(11) 
OPEN(UNIT=11,FILE=FECDF.OUT,STATUS=UNKNOWN ) 
WRITE(11,1010) NUMF,N,EB 
DO 70 I=1,NUMF 
WRITE(11,110) I,DF(D,FCD 
CONTINUE 
CLOSE(11) 
OPEN(UNIT=12,FILE=’SECDF.OUT STATUS=UNKNOWN’) 
DO 80 I=1,NIDP1 
WRITE(12,110) I,D(1),FS() 
CONTINUE 


PSPSPS SSLSE SSS SL SESS SSS ESL S SS ST SFT ST TFT SS TPS ST FT TTT FT FT ST TFT TT FT TT ST TTT TS SS 


* FORMATS 
*x 


2K KK OK ok ok ac ie oe oc oc oe oe ok kK 2 2K oo 2 oo 2 oo oe a a ok oe ae ae oe 2k ae oe oe kK 2 2K Kc eo oC aK cK KK KK OK RK KK 


100 
110 
1001 
1000 
1010 


FORMAT(1X,15) 
FORMAT(1X,15,2X,F8.4,2X,F8.4) 
FORMAT(1X,15,2X,F8.4,2X,F8.4,2X,E 10.4) 
FORMAT(1X,15,2X,E16.10,2X%,E 16.10) 
FORMAT(1X,15,2X,15,2X,F8.4) 


eke oko oc og ae oe a ke 2 2 2 aK ok a a oc 2 oe a a a oe 2 a a oe 2 ok ac ke oe ok ae ac oe oe ake 2K 2 ook 2K a ok ie ok oe KK 2K 2K ok Kk ok 


B. 


110 
303 


END 


MAXIMUM LIKELIHOOD ESTIMATOR 


PROGRAM MLE 

INTEGER NIA 

PARAMETER (NIA = 6000) 

REAL FE(NIA),DFE(NIA),FUB(NIA),DFUB(NIA),GL 
CHARACTER*1 Q1 


PRINT *, Do you want the upper bound or complete data’ 
PRINT *, ’set, suggest the compete data set when the’ 
PRINT *, "DATA is less the the number of filements’ 
PRINT *, TYPE "U" for upper and "C" otherwise’ 
FORMAT\A1) 

READ (*,303) Q1 

IF(Q1 .EQ. ’U’) THEN 
OPEN(UNIT=11,FILE=UBFECDF.OUT ,STATUS=’0LD’) 
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ELSEIF (Q1 .EQ. ’C’) THEN 
OPEN(UNIT=11,FILE=FECDF.OUT ,STATUS=’0OLD’) 
ELSE 
PRINT *, "TRY TYPING CAPITAL LETTERS’ 
GOTO 110 
ENDIF 
. 
READ(11,1010) M,N,EB 
1010 FORMAT\1X,15,2X,15,2X,F8.4) 
DO 10 I=1,M 
READ(11,1020) J,DFUB(I),FUB(I) 
10 CONTINUE 
1020 FORMAT\(1X,I5,2X,F8.4,2X,F8.4) 
3 KK ok 2 2 2 2 2 ok a 2 a 2k ok KK 2K 2 2 aK oi 2 2 2 a aK a i 2k 2h 2k 2k a a 2k 2 2 ke a 2 aK oe ae a 2 eo a ke ac ak a 2 oie ak 2c 2k 2k ok ok 
* CALCULATE ALPHA AND BETA 
PRINT *, "What is the GAGE LENGTH?’ 
READ *, GL 
* 
* REDUCE THE DATA TO ABOUT 50 POINTS 
888 PRINT *, "How many points for the MLE ?” 
READ *, MLE 
PRINT *, MLE 
DEL = 1.0(MLE+1.0) 
= 
et. 
* find the value of F above and below the expected rank and average 
400 FE(I) = 1.0*I(MLE+1.0) 
IF(I .GT. MLE) GOTO 450 
410 IF(FUBW) .LE. FE(I)) .AND. (FUB(J+1) .GT. FE(I))) THEN 
DFE(D= -(DFUB(J+1)-DFUB(J))*(FUB(J+1)-FE(I))AFUB(J+1)-FUBUWJ)) 
+ + DFUB(J+1) 
IJ=I+1 
GOTO 400 
ELSEIF((J+1) .LE. (M)) THEN 
J=J+1 
GOTO 410 
ENDIF 
450 PRINT *, I,J,FUB(J), FUB(J+1),FE(D 
IF(FUB(J+1) .LE. 0.0) THEN 
PRINT *, MLE to large try again’ 
GOTO 888 
ENDIF 
AK a 2a oie ee A a AK a Ke 2 a aK 2 2 a a aK Kok a 2 oie ok 2K a 2 2 2 ok 2K aK oc 2 2K ok aK oe 2 a 2 2 2K 2 2 2 2 2 2K 2K KK 
* estimate the initial alpha 
PRINT *, Estimate the intial alpha’ 
READ *, ALFAO 


A he 2 2 ie ie 0 2 2 2 2 2 2 a oe 2 aK oie 2 2 2 2c ok 2 2 oc oie oie oie 2 2 2 ok og 2 oC oc 2 2 og og 2 ic a 2 2 2 2c a 2 2k 2 ok 2 a ok 2k a ok 


* calculate 
ALFA = ALFAO 
KK=0 

320 DO 310 IT = 1,2 
SUM1=0.0 
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SUM2=0.0 
SUM3=0.0 


* Calculate the variables of alpha 


300 


310 


340 
341 


c. 


1010 


10 
1020 


DO 300 I=1,MLE 

IF(DFE(I) .GT. 0.0) THEN 

SUM1 = SUM1 + ((DFE(I))**ALFA * LOG(DFE(D)) 
SUM2 = SUM2 + ((DFE(I))**ALFA) 

SUM3 = SUM3 + LOG(DFE(I)) 

ENDIF 

CONTINUE 

IF(IT .EQ. 1) THEN 

ALFA1=1.0(,.SUMI1/SUM2 - (SUM3/MLE)) 

ALFA = ALFAI 

ELSEIF(UIT .EQ. 2) THEN 
ALFA2=1.0(SUMI1/SUM2 - (SUM3/MLE)) 

ENDIF 

CONTINUE 

ALFAOP = ALFAO -(ALFAO - ALFA1)**2)(ALFAO - 2*ALFAI1 


+ + ALFA2) 


TEST = ABS(ALFAOP-ALFAO) 
IF(TEST .GT. .0001) THEN 
IF(KK .GT. 1000) GOTO 341 
KK = KK + 1 

ALFAO = ALFAOP 


PRINT *, ’ALFAO =’,ALFAO 


GOTO 320 
ELSE 
ALPHA = ALFAOP 
ENDIF 
SUM = 0.0 
DO 340 I = 1,MLE 
SUM = SUM + (DFE(I))**ALPHA 
CONTINUE 
BETA = (SUM/(MLE))**(1/ALPHA) 
PRINT *, "ALPHA = ’,ALPHA,’ BETA = ’,BETA/GL, mm/mm’ 
END 


WEIBULL LINEAR CDF 


PROGRAM LIN 

INTEGER NIA,AREA 

PARAMETER (NIA = 6000) 

REAL FE(NIA), DFE(NIA),FUB(NIA),DFUB(NIA),FS(NIA),D(NIA) 


OPEN(UNIT=11,FILE=UBFECDF.OUT,STATUS=0LD’) 
READ(11,1010) M,N,EB 


FORMAT(1X,15,2X,15,2X,F8.4) 


DO 10 I=1,M 
REAIX(11,1020) J, DFUB(I), FUB(I) 


CONTINUE 
FORMAT(1X,15,2X,F8.4,2X,F8.4) 


KE 2 2K ie ae oe te 2 2 2 2 IK KC 2 2 oC 2 aie oe ae ee 2 2 KC 2 2 2 oe 2 ok kc 2 of 2 2 2 2 2 oe Kc 2 aie oe ak Cae KK 2K aK aK ok 2 2K OK ok 
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* REDUCE THE DATA TO ABOUT 50 POINTS, USING THE EXPECTED 


RANK. 
PRINT *, “HOW MANY POINTS TO FILTER USING THE EXPECTED 


READ *, NUM 
PRINT *, NUM 
DEL = 1.0({NUM+1.0) 
a | 
i 
* find the value of F above and below the expected rank and average 
400 FE(I) = 1.0*IANUM+1.0) 
IF(I .GT. NUM) GOTO 450 
410 IF((FUB(J) .LE. FE(I)) .AND. (FUB(J+1) .GT. FE(I))) THEN 
DFE(I= -(DFUB(J+1)-DFUB(J))*(FUB(J+1)-FE(I))(FUB(J+1)-FUB(J)) 
+ + DFUBUJ+1) 
J=I+1 
GOTO 400 
ELSEIF((J+1) .LT. (M+1)) THEN 
J=J+1 
GOTO 410 
ENDIF 
450 PRINT *, I,J 


ot 
KK KK KKK KKK KKK KKK KKK KKK KKK KKK KKK KK KKK KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK SK 


* CALCULATE ALPHA AND BETA 
* 


* FIND F*=LOG(-LOG(1-F)) VS LOG(DELTA/BETA) 
DO 20 I=1,M 
IF(FUB(I) .LE. 0.0) FUB(I) = 0.0001 
FS(I) = LOG(-LOG(1-FUB(1))) 
IF(DFUB(I) .LE. 0.0) DFUB(I) = 0.0001 
IXI) = LOG(DFUB(I)) 
20 CONTINUE 
DO 70 I=1,M 
510 =IF(FS(I) .LE. 0.63212) .AND. (FS(I+1) .GT. 0.63212)) THEN 
BETA = -(D(I+1)-D(1))*(FS(1+1)-0.63212)(FS(1+1)-FS(1)) 


+ + [TXI+1) 
GOTO 71 
ENDIF 

70 CONTINUE 
71 CONTINUE 


BETA = EXP(BETA) 
PRINT *,’ BETA = ’,BETA 


DO 80 I=1,M 
IF(DFUB(I) .LE. 0.0) DFUB(I) = 0.0001 
D{I) = LOG(DFUB(I)/BETA) 
80 CONTINUE 
BETA = EB*BETA/N 
Ltt L ot Sete te et eet ee Te ett T eT etre Tere Tree Teer eer ere LS eee eee Tes 
OPEN(UNIT=11,FILE=’FSTAR.OUT,STATUS=’ UNKNOWN’) 
DO 60 I=1,M 
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60 


WRITE(11,1020) I,D(1), FS) 
CONTINUE 
END 
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APPENDIX F 


SIMULATION FORTRAN CODE 
PROGRAM SIMUL 


- SIMUL IS A STATISTICAL SIMULATION ROUTINE USED 
* FOR RELIABILITY 
- ANALYSIS OF COMPOSITE FILEMENTS. 
UTILITIES FILES REQUIRED ARE: 
SEED - file containing the integer seed for random 
number generation 
INPUT - file containing the input parameters * 


OUTPUT FILES INCLUDE: 
EXPER.OUT - the experimental simulation of discrete 
data. 
SEED - outputs a new seed. 


* * &* &€& & HH HN HR OH OF 


* INPUT FILE FOR BUNDLE2.FOR PROGRAM __ siline 1 
* XXXXXXXKXKXXKXXXKKKKK 


12345678 ‘line 2 
* E (kg) .90 cline 3 
* ROWS = M 100 ‘line 4 
* COLS = N 100 ‘line 5 
*MGL (mm) 250.0 sline 6 
* ALPHA 4.0 ‘line 7 
* BETA 16.0 ‘line 8 
* MODEL NORMAL ‘line 9 
* SKEW/ALPHA 3.5 ‘line 10 
* MAX SLACK (% _ ) 10.0 line 11 % OF MAX 
‘3 DISPLACEMENT 
* SAMPLE RATE 20.0 ‘line 12 pts/s 
* CROSS-HEAD SPEED 1.8 ‘line 13 mm/min 
* NOISE +/- kg .001 sine 14 tolerance 
* TOLD (mm) .00254 ‘line 15 tolerance of displ. 
. output 


* 


* INPUT PARAMETERS ARE: 

: 1) THE LENGTH OF THE TEST BUNDLE (MGL). 

i 2) THE NUMBER OF FILEMENTS PER BUNDLE (N). 

me 3) THE NUMBER OF SEGMENTS PER FILEMENT (M). 

: 4) THE SEGMENT LENGTH IS DEFINED BY MGL=M*LS 
* 
* 


THIS PROGRAM CREATES: 
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1) THE RANDOM NUMBERS 

2) CONVERTS TO WIEBULL BUNDLE 

3) FIND THE WEAKEST IN THE ENTIRE BUNDLE 
FOR THE LOWER BOUND. 

4) FIND THE WEAKEST IN EACH FILEMENT 

5) ORDER THE WEAKEST TO STRONGEST FILEMENTS 

6) CHANGE THE STRENGTHS TO BREAKING LOADS AND 
DISPLACEMENTS 


WK KK KKK KK KKK KKK KEKE KKK KKK KEK KK KEKE KKK KEKE KKK KKK KEKE RKKKKEKEKK 


* Begin Program: 
* 


* & He HH HH HH H 


- SSeS SS SSS SS SSS SELES SSS SSS SSS LS SS SSS STS FS SST SSS SSS SSS SFT SSS ST SS FS 3 
* 


INTEGER ROWS,COLS,MAX,LIM,LIM1,N,PTS,COUNT,LIM2 

PARAMETER(LIM=6000,LIM1=1000, LIM2=1001) 

REAL BUN(LIM1),RAND1(LIM1),RAND2(LIM1),DISPIXLIM1), 
+SLACK(LIM1),P(LIM1),ALPHA,BETA,MGL,BRKLDM, 
+DISPLM,MAXX MAXY,MEAN,PU(LIM1), 
+PI{LIM1),SUM,MAXSLK,TEMP(LIM2),PS(LIM2),D(LIM), 
+DISPLA(LIM1),DEL,PMS,LOAD(LIM),NSE 

DIMENSION NO(LIM1),NB(LIM1),NDB(LIM1) 

DOUBLE PRECISION SEED 

CHARACTER*1 MODEL*7,SONS,Q,Q2,Q3,Q4 

COMMON MAX 

COMMON /SML/M,N 


* 
KEKKKKKKKKKKK KKK KKK KKK KKKK KK KKK KKKKK KKK KK KKKKKKKKKKKKKEKKKKKKKKEKEKKE 
* 


* FILE 11 = PARAM DATA FOR INPUT 
* FILE 12 = SEED DATA FOR SEEDING 


« 
KK KK KKK KKKKKKK KKK KKK KKKK KKK KK KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKEKEKKK 
~ 


OPEN(UNIT=11,FILE=INPUT,STATUS=’0OLD’) 
OPEN(UNIT=12,FILE=’SEED’,STATUS=’OLD’) 
* 


HE EE HH fe FH 2 oe 2 A EK A EK EK KK KKK KKK KK KK 


ig READ IN THE INPUT DATA 
* 


KKK KK KKK KK KKK KKK KKK KK KKK KKK KKK KKK KKK KKK KKK KKKKKKKKKKKKKEKKEK K 
x 


READ\(11,1530) 
REAIX(11,1500) E 

PRINT *, MODULUS=2’, E 
REAIX(11,1520) M 

PRINT *, M = ’,M 
ROWS=M 

REAIDX(11,1520) N 
COLS=N 

PRINT *, N = ’,N 
REAIX(11,1510) MGL 
PRINT*, ’Mean Gage Length = ’,MGL 
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MAX=M*N 
READ(11,1510) ALPHA 
PRINT *, “ALPHA = ’,ALPHA 
READ(11,1510) BETA 
PRINT *, "BETA =’,BETA 


* 
* * * * ae * 


* Determine which model to use for slack. 
* * * * * © 


* 


REAIDX(11,1515) MODEL 
PRINT *, "MODEL=’,MODEL, APPROXIMATION’ 
IF (MODEL .EQ. "NORMAL’) THEN 
READ(11,1500) STDDEV 
PRINT *, "ALPHA = ’,STDDEV 
READ(11,1510) MEAN 
PRINT *, "Maximum slack as percent of expected’ 
PRINT *, ’maximum displacement (3% gage length).=’,;MEAN 
MEAN = (MEAN/100)*(0.03)*MGL/2 
PRINT *, Mean=’,MEAN,’mm’ 
PRINT *, "Maximum slack =’,MEAN*2,mm’ 
ELSE 
READ(11,1500) RANGE1 
REAIX(11,1500) RANGE2 
RANGE 1=(RANGE1/100)*(.03)*MGL 
RANGE2=(RANGE2/100)*(.03)*MGL 
PRINT *, "RANGE 1 =’,RANGE1,’ mm’ 
PRINT *, "RANGE 2 =’,RANGE2,’ mm’ 
ENDIF 
x 


* ae sx Bs * %e 
* Read in the sampling rate and cross-head speed. 
* aR K * * % 


Ae AR A ee 6 2 ok a ie A 2 2 2 2 oR OK aK 2 aK 2 2 2 2K oR 2 aK oR oo iC aK 2 aK a 2 2 oi oc 2 2 2 aK 2 2K aK a Co 2K aK 2K 2K ok OK A OK ok Ko ok 2k 


* Ensure the displacement interval governed by the cross-head and 
* sampling rate is greater than the tolerance level of the output 

* displacement data. 

* 


3K 2 2 2 2 2K 2K 2 i 2K 2 2 og 2 2 2 2 og aK ie oo 2 ok 2 2 oie 2 ok 2 2 2 2 2 a 2 2 oie oie ok oie 2 2 kc 2 a 2 2 2g 2 ok kc 2 2g 2 ok 2 2k ok 2k aK ok 2k ok ok kk 
* 


READ(11,1510) SR 

PRINT *, Sampling rate =’,SR 

REAIDX(11,1510) XSP 

PRINT *, ’Cross-head speed =’,XSP 
* 


£ * * xk x Xk 
* How much noise. 
7 * aK aK * ae 
* 

READ(11,1510) NSE 

PRINT *, Range of noise =’,NSE 

REAIX(11,1510) TOLD 
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PRINT *, Tolerance of output =’,TOLD 
x * % * * * 
* Calculate the number of total points recorded in the as 
* 


* experiment. 
* 3K 


‘é 


444 RR = MGL*(0.03) 
PTS = RR*SR/(XSP/60) 
DELTAD = RR/PTS 
IF(TOLD .GT. DELTAD) THEN 
PRINT *, The displacement interval is less than the output’ 
PRINT *, ’tolerance, please input a new cross-head and/or’ 
PRINT *, ’sampling rate.’ 
PRINT *, Current Cross-head = ’,XSP 
PRINT *, Input new XSP’ 
READ *, XSP 
PRINT *, Current sampling rate = ’,SR 
PRINT *, Input new sampling rate’ 
READ *, SR 
GOTO 444 
ENDIF 


* ae oK * * 


* * * * * * 


* Determine the type of simulation. 
3K * * * * 


PRINT *, "SLACK OR NO SLACK ? (Y/NY 
READ (*,1600) SONS 
PRINT *, "NOISE OR NO NOISE ? (Y/NY 
READ (*,1600) Q2 
PRINT *, ’Change the SEED ? 
READ (*,1600) Q 
PRINT *, ’?COMPLAINCE (Y/N)? 
READ (*,1600) Q3 
IF((Q3 .EQ. Y) .OR. (Q3 .EQ. ’y’)) THEN 
PRINT *, WHAT IS THE COMPLIANCE NOISE LEVEL? 
READ *, CTOL 
ENDIF 
PRINT *, Do you want English units (lbs,inches) on output? 
+(Y/NY 
PRINT *, ’If N the output is in SI units (kg,mm)’ 
READ (*,1600) Q4 


CLOSE(11) 
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AE KK AE KR A KK A OK OK EK EK OK KO KK KK OK OK KKK KKK KK KK KE 


* BEGIN THE ALGORITHM ts 
AK 2K 2 2 2 2K A KK 2K 2K 2 2 EK KK OK EK KK RK KK OK OE KE ER KK OK OK OK KE KK OK OK KK OK OK EK OK KOK OK KK OK 

* 

td * * 3 * * * * 


* Read the seed. (seed for the uniform random number generator). 
* and generate the random numbers 
+ * * * oe 


* * * * 
* 

READ(12,1000) SEED 

REWIND(12) 

DO 10 I=1,N 


CALL RANDXSEED,M,BUN) 

CALL SMALKROWS,BUN,RAND1(I)) 
10 CONTINUE 
* 


** * ae * * * * 2K 
* Check if slack is requested in simulation. 

OK * * * * aK 2K * 
* 


IF (SONS .EQ. ’Y’) THEN 
CALL RAND(SEED,N,RAND2) 


ENDIF 
* 
* * ae * * 5 x * 
* Initialize the order arrays. 
* K * * * * * oR a 
* 
DO 25 I=1,N 
SLACK(I)=0.0 
NO(D=1 
NB(I)=] 
NDB(I)=1 
25 CONTINUE 
* 
KK * * * 9 » * * 


* Check if seed is to be changed, if so write the new seed. 
* * * 5 ae 5 * 


* 
IF ((Q .EQ. Y) .OR. (Q .EQ. ’y)) THEN 
WRITE(12,1000) SEED 
REWIND (12) 
CLOSE(12) 
ENDIF 


Ae Ae A A A A AE OE IE KK KK KK RRR RRR ROKK 


* Generate the load and slack distributions. 
* 


EK EK ER KR oR KK RK KR KKK KR oR KK RRR KKK RK ROR KK KK RK aK oR oo oo ok oo 2 ok ok ok ok ok ok ok 
* 
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CALL WEIBUL(ALPHA,BETA,N,RANDI,P) 


IF (SONS .EQ. ’Y .OR. SONS .EQ. ’y’) THEN 
IF (MODEL .EQ. NORMAL’) THEN 
CALL NORMAL(STDDEV,MEAN,N,RAND2,SLACK,MAXSLK) 
ELSE 
CALL UNIFM(RANGE1,RANGE2,N,RAND2,SLACK,MAXSLK) 
ENDIF 
SLACK(1)=0.0 

ENDIF 

* 


AE AR A 2 oo KK RK eo oo 2 oi oo oo og eK oo Ko Kok ok ok ok ok ok ok ok ok ok ok ok 
* Find the displacement for each failure load and slack. 
x 


ARO RR Ko KK KK KR KK RK KK KK KK KK KKK KK KK OK KK KK KKK Kok Kok Ok 
* 


DO 20 I=1,N 
DISPLAT) = (P(D/E)*(SLACK(I) + MGL) 
20 CONTINUE 
=iN| 
* * * * * 


Order the "DISPL" array, and carry "No" 
* aK ak *x x 


* + + 


CALL SHELL(N,DISPL,NO) 


% 


* * * * * 


* reorder the remaining arrays, the others are 
keyed to "DISPL" 
* x 


%* + 


* * * 


CALL SWITCH(N,NO,RAND1) 
CALL SWITCH(N,NO,RAND2) 
CALL SWITCH(N,NO,SLACK) 
CALL SWITCH(N,NO,P) 


* 

AE OR A EE KR KK KK RK KKK KKK KKK KKK KK KK KKK 
* Calculate the apparent displacement. 

ARK 2 A i EK ie i i 2 a i i 2 OK KK KK RK RK KK KK KKK KKK KKK KK KK KKK KKK 
* 


DO 40 I=1,N 

DISPLA(I) = DISPL{I) + SLACK() 
40 CONTINUE 
* 


KK * * * * * 


* Reorder the arrays keyed to "DISPLA" 
* * aE x x 


ye 


CALL SHELL(N,DISPLA,NB) 
CALL SWITCH(N,NB,RAND1) 
CALL SWITCH(N,NB,RAND2) 
CALL SWITCH(N,NB,SLACK) 
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CALL SWITCH(N,NB,DISPL) 
CALL SWITCH(N,NB,P) 
CALL SWCH2(N,NB,NO) 

* 


* * * * * * 


* Define the ordered slack as "TEMP" 
* ok a + ak * 


* 


DO 45 I=1,N 
TEMP(I)=SLACK(T) 
45 CONTINUE 
IF((SONS .EQ. ’Y) .OR. (SONS .EQ. ’y’)) THEN 
CALL SHELLAN,TEMP,NDB) 


* 
Po PES TSS SLE SSS SSS SLT LSS SSS TT SL ST TT SST TS TST TSS TELE LLL SSS FS FL FS FSF FS +S tS FF | 
+ a 


* Calculate the maximum slack and intercept point. - 
* Calculate the upper and lower Pi bundle load at each displa. * 
a 


* 
ea * 

Re 2 2 2 OR 2 ie 2 2 2 2 ie ie ie 2 2 i eK ee 2 2 2 2 kk ie 2 2 ke 2 kk ok ee ek ee 2 2 2 2k 2k ok ok ok ok ok Ok 
* Find the loads at the slack points. 

KK * ae ik * 
* 


oF 


PS(1)=0.0 

DO 70 J=2,N 

CONST=0.0 

DO 75 [=1,(J-1) 

CONST=CONST+E*(TEMP(J)-TEM P(1))(MGL+TEMP(1)) 

75 CONTINUE 

PS(J)=CONST 
70 CONTINUE 


25 2k 2 2 2 2 2 2 2 gg 2 2K 2g 2 2 2 2K 2k 2k 2 oie 2 ok ik 2k ok ie 2 ik ik 2 ok 2k 2k oe ok 2k oie og 2 oe ok oie 2 2 ig 2k 2k oe ok 2 2k ok 2c 2k ok ok ok 2 2 2K ok ok 2k ok 


* ADD COMPLIANCE TO SLACK 
* 


% 


* DC1 = mm * 


AE A Ae OR 2 2 OR 2 ig 2K 2 2 6 2 2 2 2 2 oe oo 2 eK 2 ok 2k ie og 2 OK ek ie og ok ok i oo 2 ok io ok ak ok 


IF((Q3 .EQ. ’Y’) .OR. (Q3 .EQ. ’y’)) THEN 


DC1 = .0356 
Al= 0.0512 
A2= -0.0217 
POWER = 0.1975 
A3 = 0.0064 
A4 = 0.0242 


DO 120 I = 1,N 
IF(TEMP(1) .LT. DC1) .AND. (PS(I) .GE. CTOL)) THEN 
TEMP(I) = TEMP(I) + A1*(PS(I)**POWER) + A2 
ELSEIF(PS(I) .GE. .01) THEN 
TEMP(I) = TEMP(I) + A3*PS(I) + A4 
ENDIF 

120 CONTINUE 
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ENDIF 
* 


Prt SSS TSS SSS SSS SSS SSS SSF TFS SF SET FSF S SSS STL FSS TTF TFT SST TFT TS FFF FF FSF FF SF 
* FAILURE POINTS. 

2 2 2 2 2 oo oo oo oo ie 2 2 2 2 2K 2 2 2 2 2 2 2 2 2 2 2K KK 2 EK ok 2 OK KO ok ok ok ok KK 
* 


DO 80 J=1,N 
CONST=0.0 
DO 85 I=J,N 
CONST=CONST+E*(DISPLA(J)-SLACK(1))(MGL+SLACK(I)) 
85 CONTINUE 
PU(J)=CONST 
80 CONTINUE 


2 2 2 EK 2 ok oo oo ok 2 ok oo 2 2 oo 2 ok oe 2 oo oo Ko oo oR Kok ok 


* ADD COMPLIANCE TO FAILURE REGION 
* 


* 


*2DC y=) min - 
- ES SS SSS SSS SES SSS SSS SFT SSS PSPS SSS SSS SPS SS SLL SFT SSS STS STS FT FSS SS SS SSS SS SS SS 
IF((Q3 .EQ. Y) .OR. (Q3 .EQ. ’y’)) THEN 
DO 121 1 = 1,N 
IF((DISPLA(I) .LT. DC1) .AND. (PU(I) .GE. CTOL)) THEN 
DISPLA(I) = DISPLA(I) + A1*(PU(I)** POWER) + A2 
ELSEIF(PU(I) .GE. .01) THEN 
DISPLA(I) = DISPLA(I) + A3*PU(I) + A4 
ENDIF 
121 CONTINUE 
ENDIF 


* 
2K KK 2 2 2 oo 2 2 2 2 2 2 2 2 KE EK RO OK KKK KK EK 


* ADD DISCRETE DATA * 
* TO THE * 

* SLACK REGION. * 

* 


* 
-S SPSS SSS SSL SS SSS SSS SSSSS SSS SSS SESS SS SSS SSS SS SPSS LSS SSS LSS SSS SS SS SS Ft | 
* 
J=1 
ADD = 0.0 
CONST=0.0 
D(J)=DELTAD 
DO 90 I=1,(N-1) 
IF (DJ) .LT. TEMP(I+1)) THEN 
SLOPE = (PS(I+1)}-PS(1))KTEMP(I+1)-TEMP(])) 
91 LOADXJ)=(IX J)-TEMP(I))*SLOPE + ADD 
CONST=DJ) 
J=Js+1 
D(J)=CONST+DELTAD 
IF (XJ) .LT. TEMP(I+1)) GOTO 91 
ENDIF 
ADD = ADD + PS(I+1)-PS(1) 
90 CONTINUE 
JBIGL = J-1 
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OK oe kK a oc 2h 2 2k 2 2 2 2 oi oie oe a oi 2k 2k ik 2c 2 2 2 2c 2c 2c 2 2 aK aK oe cc EC EK ok 2 OK ok 2 oo a a aK oi 2 2 2 2c 2K KK KK ok oR OK 


* ADD DISCRETE DATA * 
* TO THE * 

* LINEAR REGION * 
rs 


x 
3K 2K ok 2k 2k 2 oe i 2 2 2 2 Ko IK 2K 2 2 i oo eK 2 2 OK 2 2K 2K 2 2K 2 2K 9 2K 2K 2K 9 oi 2 2K 2k 2 9 9 2k 2k 9k 2K oi oko 2k 2K 2 OK 2k KK OK ok ok ok kK 
Ps 


CONST=0.0 

SLOPE=(PU(1)-PS(N)(DISPLA(1)-TEMP(N)) 
100 LOAD(J)=SLOPE*(D(J)-TEMP(N)) + ADD 

CONST= DJ) 

J=J+1 

D(J)=CONST + DELTAD 

IF (D(J) .LT. DISPLA(1)) GOTO 100 

JENDL = J-1 


* 
* * * * * * 


* Find the intercept point. 
£ * ae ae 


* 
DINT= DISPLA(1) - PU(1)/SLOPE 

* 

2 2K oie ok 2 2 2 2k 2 oi ok 2 2 2K ok ok 2K 2 2K 2k 2 2 aK 2K 2K 2 aK 2K 2K 2K 2k 9K oR 2k 2K oo 2 Kok 2 ok 2 oi 2 2K ok KK ok ok 2K ok ok ok ok oR Ko KOK 


* ADD DISCRETE DATA - 
* TO THE * 

* FAILURE REGION. 
* 


* 
OE 2 2 2 a 2 2 i 2 aK 2 2 2K i 2 2 2 2 aK aK 2 2 OK oi 2c 2c 2 2 2k 2 2c oe 2 2k 2k 2k 2 a oo aK 2k 2k 2k 2k 2k 2k 2k 2 2 2 2k ok ok ok ok ok 2k 


DO 110 I=2,N 

IF (D(J) .LT. DISPLA(I)) THEN 
SLOPE = PU(I)ADISPLA(I)-DINT) 

rs LOAD(J)=SLOPE*(D(J)-DINT) 

CONST = D(J) 
J=J+1 
D(iJ) = CONST + DELTAD 
IF (DiJ) .LE. DISPLA(I)) GOTO 115 


ENDIF 
110 CONTINUE 
LOAD(J)=0.0 
700 PRINT *, "CONTINOUS DATA COMPLETE ’ 
COUNT = J 
PRINT *, "COUNT=’,J 
* 
* * * * * * 
* Find the largest value for plotting 
* * * * * * 
* 
CALL LARGE(K,PU,BRKLDM) 
« oe * * * * 
* Where brkldm and sldm are the maximum for that vector. 
* * * * * * 
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DISPLM=DISPLA(N) 


Re RE 2 2 2 2 a a 2c 2k ok ok ic 2c 2k oi aK oe ke ke ke 2 ok akc ake ok oe ake 2 ec kc oe oc kc oe oe oc kc ok ok 2 2 2K oe ok 2c oc ok oo a ook ok ok 


* CREATE THE NOISE 
* 


26 2K 2K oe oc 2 2k oc ke 2k ok ak kc 2c 2k ok ok 2k ek ok aK oe oc oe 2k ake ke of 2 2 oe ok oe ke 2 oe ake 2 ok oe 2K oe oe kc ke 2K ke ok 2K ok oie oe oe kc 2c ok aK oe ke ok ok ok ok 


* 
IF ((Q2 .EQ. 'Y’) .OR. (Q2 .EQ. ’y’)) THEN 
CALL NOISE(SEED,COUNT,LOAD,NSE) 
ENDIF 
* 
PSS SSFP LSE LSS SSS SLES ESS SSE SES SSS RST SE SESE SESE SSESLS SLE TSES ELE LESS SSS TF + SF | 
* PRINT THE DATA *: 
* * 


2H 2 2 2 oie oe ke 2 2 oe 2c 2k ok ak ke 2K ke 2k 2k Kk oe ake 2k oe 2K eke 2k ak ak kek ak a 2c 2g ok age ak a a 2c oe ak ok a oe a ake ok ok ac oc ok ok oe ok ok ok ok 


OPEN(UNIT=13,FILE=’"EXPER.OUT’ ,STATUS=’UNKNOWN’) 
IF((Q4 .EQ. Y’) .OR. (Q4 .EQ. ’y)) THEN 

CONST1=0.03937 

CONST2=2.2046 


ELSE 
CONST1 = 1.0 
CONST2 = 1.0 
ENDIF 


DO 105 Il=1,COUNT 

WRITE(13,111) I,.D()*CONST1),(LOAD(I)*CONST2) 
105 CONTINUE 

WRITE(13,116) ALPHA,BETA,N 

WRITE(13,111) COUNT,DISPLM,BRKLDM 


CLOSE(13) 
TrTtttt ttt et te TTT TTT TTT TTT TT Te Te TT TTT TTT TT TTS tet ee eee Se eee ee 
* FORMATS is 
* _ 


2 2k 2k 2k a 2k 2 a og 2 2c 2 oe 2 a ke ok oe ke 2k ok oe 2c ke ok a 2 aK oe ke 2 ok 2c 2 ok oe oe ke ake 6 2k Kg 2c 2 ak i 2 ok 2 2 ok ok ok 2 ok ok ok OK ok ok ok 


111 FORMAT(1X,15,2X%,F8.4,2X ,F8.4) 

116 FORMAT (1X, PARAMETERSALPHA2]’,F4.1,2X, BETA’, 
+F4.1,/,1X% 
+, N=’,14,/,1X,’N’,5X, DELTA’,10X, LOAD’) 

1000 FORMAT\(1X,F15.1) 

1500 FORMAT(20X,E16.10) 

1510 FORMAT(20X,F8.4) 

1515 FORMAT(20X,A6) 

1520 FORMAT(20X,[I5) 

1530 FORMAT(1X,/) 

1600 FORMAT\A1) 

* 


STOP 
END 


SUBROUTINE NOISE(SEED,COUNT,LOAD,NSE) 
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x 


gee ee 


% pd 


INTEGER COUNT 

DOUBLE PRECISION SEED 

REAL LOAD(COUNT),NSE,RNDM 

DO 10 J=1,COUNT 

CALL RANDN(SEED,RNDM) 

LOAD(I) = LOAD(I) + (RNDM*NSE - NSE/2) 
CONTINUE 

END 


INTEGER N,ORDER(N),LIM1 

PARAMETER(LIM=5000) 

REAL ARRAY(N),TEMP(LIM) 

DO 1 I=1,N 
TEMP(1I)=ARRAY(T) 
CONTINUE 


DO 2 I=1,N 
K=ORDER(]) 
ARRAY(I)=TEMP(K) 

CONTINUE 

END 


INTEGER N,LIM1 
PARAMETER(LIM1=5000) 
INTEGER ARRAY(N),TEMP(LIM1),ORDER(N) 
DO 1 I=1,N 

TEMP(T)=ARRAY(T) 

CONTINUE 


DO 2 I=1,N 
ARRAY(I)=TEMP(ORDER(])) 

CONTINUE 

END 


95 


* * * 


THIS ROUTINE TAKES A MATRIX IN VECTOR FORM AND 
FINDS THE SMALLEST 
VALUE IN EACH ROW AND PUTS IT INTO A VECTOR 
INTEGER TOTAL,ROWS,COLS 
REAL VECTOR(TOTAL),LOW 
COMMON /SML/ ROWS,COLS 


INTEGER MAX 

REAL VECTOR(MAX), MAXVAL 
MAXVAL=VECTOR(1) 

DO 10 J=2,MAX 

IF (VECTORUJ) .GT. MAXVAL) THEN 
MAXVAL = VECTOR(J) 

ENDIF 

CONTINUE 

END 


* LOOK AT ONLY ROWS 


* 


1 


LOW=VECTOR(1) 


DO 10 I=1,ROWS 
IF (VECTOR() .LT. LOW) LOW=VECTOR() 


CONTINUE 


END 


MAXNUM IS THE NUMBER IN THE VECTOR TO BE 
SORTED 


INTEGER MAXNUM,OFFSET,SW,NO,LIMIT 
DIMENSION NO(MAXNUM) 

REAL ARRAY(MAXNUM),TEMP,TEMP2 
LOGICAL SWITCH 


PRINT *, "SUB SHELL’ OFFSET=MAXNUM/2 
IF (OFFSET .GT. 0) THEN 


LIMIT=MAXNUM-OFFSET 


D 


* 


SWITCH= .FALSE. 


DO 15 J=1,LIMIT 
IF (ARRAY(J) .GE. ARRAY(J+OFFSET)) THEN 

TEMP = ARRAY(J) 

TEMP2 = NOG) 
ARRAY(J)=ARRAY(J+OFFSET) 
ARRAY(J+OFFSET)=TEMP 

NO(J+OFFSET)=TEMP2 
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15 


SWITCH = .TRUE. 
SW=J 
ENDIF 
CONTINUE 
LIMIT = SW -OFFSET 
IF(SWITCH) GOTO 5 
OFFSET=OFFSET/2 
GOTO 1 
ENDIF 


END 


Random number generator - Uniformly distributed in (0,1) 
SEED in (1,247483647) 

DOUBLE PRECISION RANDM,ASEED 

INTEGER TOTAL 

REAL RNDM(TOTAL) 


DO 20 I=1,TOTAL 

ASEED=ASEED+1.0 

RANDM = DMOD(16807.0D0*ASEED,2147483647.0D0) 
ASEED=RANDM 

RNDM(I) = SNGLK(RANDM*4.6566128752458D-10)-1.0E-07 
CONTINUE 

ASEED=ANINT(ASEED) 


Random number generator - Uniformly distnbuted in (0,1) 
SEED in (1,247483647) 

DOUBLE PRECISION RANDM,ASEED 

INTEGER TOTAL 

REAL RNDM 

ASEED=ASEED+1.0 

RANDM = DMOIDX(16807.0D0*ASEED,2147483647.0D0) 

ASEED=RANDM 

RNDM = SNGL(RANDM*4.6566128752458D-10)-1.0E-07 

ASEED=ANINT(ASEED) 

END 


SUBROUTINE NORMAL(ALPHA,BETA,N,RAND2, 


+SLACK,MAXSLK) 


PERFORMS A SIMULATION OF EXPERIMENTS USING THE * 
WEIBULL MODEL 


INTEGER N 
REAL ALPHA,BETA,MEAN,MAXSLK,RAND2(N), 
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+SLACK(N),RNDM 
MAXSLK = 0.0 
PRINT *, "SUB NORMAL’ 
*.. CALCULATE THE STRENGTH OF THE WEIBULL FUNCTION, 
DO 10 I=1,N 
RNDM = RAND2(I) 
SLAC K(I) 
+EXP((LOG(-LOG(1- RNDM)}+ ALPHA*LOG(BETA)VALPHA) 
IF (SLACK(I) .GT. MAXSLK) THEN 
MAXSLK = SLACK(T) 
ENDIF 
10 CONTINUE 
END 
" END 
SUBROUTINE 
+ UNIFM(RANGE1,RANGE2,N,RAND2,SLACK,MAXSLK) 
INTEGER N 
REAL RANGE,MEAN,RAND2(N),SLACK(N),MAXSLK 


CONST=RANGE2-RANGE1 
MAXSLK=0.0 
DO10 JEEN 
SLACK(I)D=RAND2(1)*CONST + RANGE] 
IF (MAXSLK .LT. SLACK(I)) THEN 
MAXSLK=SLACK(T) 
ENDIF 

10 CONTINUE 
END 


PERFORMS A SIMULATION OF EXPERIMENTS USING THE 
WEIBULL MODEL 
INITIALIZE 
REAL ALPHA,BETA 
INTEGER TOTAL,MAX 
REAL P(TOTAL),BUNDLE(TOTAL),RNDM 
PRINT *, "SUB WEIBUL’ 
*.. CALCULATE THE STRENGTH OF THE WEIBULL FUNCTION, 
DO 10 I=1, TOTAL 
RNDM = BUNDLE(]) 
rp 
+ EXP((LOG(-LOG(1-RNDM))+ALPHA*LOG(BETA))/ALPHA) 
10 CONTINUE 
END 


* + * 
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